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Applied Differential Equations 


A Course for Science and Engineering 


Organization. [Each chapter of the text is organized in sections that 
represent one or two classroom lectures of 50 minutes each. The outside 
work for these divisions requires one to six hours, depending upon the depth 
of study. 


Each section within a chapter consists of three distinct parts. The divi- 
sions represent the lecture, examples and technical details. Generally, 
proofs of theorems or long justifications of formulas are delayed until after 
the examples. The lectures contain only the briefest examples, figures and 
illustrations. 


A key to a successful course is a weekly session dedicated to review, drill, 
answers, solutions, exposition and exam preparation. While group meetings 
are important, individual effort is required to flesh out the details and to 
learn the subject in depth. The textbook design supports targeted self-study 
through its examples and exercises. 


There is a defense for this style of presentation, matched equally by a long 
list of criticisms. The defense is that this style represents how material is 
presented in classroom lectures, and how the topics are studied in the private 
life of a student. Certainly, students don’t read everything in the textbook, 
and this style addresses the issue of what to skip and what to read in detail. 
The criticisms include a departure from standard textbooks, which intermix 
theory and examples and proofs. Additional criticisms include a general need 
to flip pages to look up details. 


Prerequisites. Beginning sections of chapters require college algebra, 
basic high school geometry, coordinate geometry, elementary trigonometry, 
differential calculus and integral calculus. Several variable calculus and linear 
algebra are assumed for certain advanced topics. Instructors are the best 
judges of what to include and what to skip, concerning advanced topics in 
this textbook. 


Survey course. A complete survey course in differential equations for 
engineering and science can be constructed from the lectures and examples, 
by skipping the technical details supplied in the text. Interested students 
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can read the details to obtain a deeper introduction to the subject. Such 
survey courses will necessarily contact more chapters and trade the depth of 
a conventional course for a faster pace, easier topics, and more variety. 


Conventional Course. Differential equations courses at the under- 
graduate level will present some or all of the technical details in class, as 
part of the lecture. Deeper study with technical details is warranted for spe- 
cialties like physics and electrical engineering. Hybrid courses that combine 
the conventional course and the engineering course can be realized. 


To the Student. Expertise in the background topics is expected of 
students only after review and continued use in the course, especially by 
writing solutions to exercises. Instructors are advised that an exercise list 
and subsequent evaluation of the work is essential for successful use of the 
text. 


Matched in the text are examples, exercises and odd answers. To learn the 
subject, not only is it required to solve exercises, but to write exercises, 
which is not different from writing in a foreign language. 


Writing requires two or more drafts and a final report or presentation. En- 
gineering paper and lineless duplicator paper encourage final reports with 
adequate white space between equations. Pencil and eraser save time. Pens 
and word processors waste time. 


Contributions to legibility, organization and presentation of hand-written 
exercises were made at The University of Utah, by numerous creative en- 
gineering students, over the years 1990-2016. Their ideas produced the 
suggestions below, which were applied to the text examples and illustra- 
tions. 


1. A report is hand-written by pencil on printer paper or engineering 
paper. It starts with a problem statement followed perhaps by a final 
answer summary. Supporting material appears at the end, like a tax 
return. 


2. Mathematical notation is on the left, text on the right, often a 60% 
to 40% ratio. One equal sign per line. Justify equations left or align 
on the equal signs. Vertical white space separates equation displays. 


3. Text is left-justified on the right side. It includes explanations, refer- 
ences by keyword or page number, statements and definitions, refer- 
ences to delayed details (long calculations, graphics, answer checks). 


4. Every report has an answer check. It is usual to see back of book 
as the only detail. Proofs have no answer check. 


5. No suggestion is a rule: invent your own style by theft of good ideas 
and rejection of unsuitable ideas. 


CONTENTS iii 


Work and School. Students studying in isolation are increasing in 
number, because their jobs and family drive their university schedules. In 
spite of their forced isolation from the classroom, working students with 
families seek advice from their instructors by telephone, email and office 
visits. They make use of study groups and supplemental instruction. Course 
design in universities has evolved to recognize the shift from a predominantly 
non-working student population to its present constituency. 
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Exercises in Progress, August 2016. 
This PDF is a draft of my textbook written over the years 1999-2016. 
Please, do not distribute this PDF, because it contains many errors, 
as yet undiscovered. 


Ch 1. Already completed 74 + 52 + 50 + 46 4+ 44 + 56 = 
322. 


Ch 2. Already completed 86 + 64 + 64 + 70 + 66 + 98 + 
58 + 42 + 30 + 70 = 648. 


Ch 3. Already completed 36 + 80 + 50 + 48 + 22 = 236. 


Ch 4. Already completed 52 + 38 + 34 + 24 + 20 + 28 + 
24 + 40 + 12 = 272. 


Ch 5. Already completed 128 + 63 + 106 + 90 + 58 = 445. 


Ch 6. Already completed 87 + 36 + 0 + 68 + 30 + 63 + 
36 + 75 + 46 = 278. Need repair for missing exercises in 6.3 
and 6.9 


Ch 7. Already completed 30 + 144+ 20+ 84+ 164 51 = 
139. Add 10 to 7.4. 


Ch 8. Already completed 26 + 68 + 16+0+4+0+0+0+4 
0 = 110. Add as follows: 8.3 += 40, 8.4 += 40, 8.5 += 40, 
8.6 += 30, 8.7 += 30, 8.8 += 10. 


Ch 9. Already completed 35 + 0 + 66 = 101. Add as follows: 
9.2 += 40. Fix blanks in 9.1, 9.3 


Ch 10. Already completed 0 + 0+ 0+0=0. Add 30 to 
each of 10.1, 10.2, 10.3 and 20 to 10.4. 


Ch 11. Already completed 0+ 0+0+4+544+0+94+0+4 
0 + 0 = 148. Add as follows: 11.1 += 30, 11.2 += 40, 11.3 
+= 40, 11.4 has blanks, 11.5 += 30, 11.7 += 20, 11.8 += 
20, 11.9 += 30. 


Ch 12. Already completed 0 +0+0+0+0+10+4+4+4 
26 = 40. Add as follows: 12.1 += 10, 12.2 += 40, 11.3 += 
30, 11.4 += 30, 11.5 += 30, 11.6 += 20, 11.7 += 20, 11.8 
has blanks. 


Appendix: Already completed 54 + 38 + 28 + 24 = 144. 


About 2883 problems are already prepared. More to come, 
about 740. 


CONTENTS 


Indexing 
Did ch 1,2,3,5 
To do: other chapters. Record here when finished. 


Chapter 11 


Systems of Differential 
Equations 


Contents 
11.1 Examples of Systems .............. 740 
11.2 Basic First-order System Methods ...... 765 
11.3 Structure of Linear Systems .......... 773 
11.4 Matrix Exponential ............... 781 
11.5 The Eigenanalysis Method ........... 789 
11.6 Jordan Form and Ejigenanalysis ........ 797 
11.7 Nonhomogeneous Linear Systems....... 811 
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Linear systems. 


tions of the form 


(1) 


where’ = d/dt. Given are the functions a;;(t) and f;(t) on some interval 


/ —_ 

Vy = Al2 sot AiniIn T+ fi, 
! 

Tg = Qe aoe A2ntn fa, 
/ = 

Lm = Am1i%1 ois Amntn 1 Fis 


a<t<b. The unknowns are the functions x(t), ..., ¢n(t). 


The system i 


s called homogeneous if all f; = 0, otherwise it is called 


non-homogeneous. 


Matrix Notation for Systems. 
linear equations (1) is written as the equivalent vector-matrix system 


x’ = A(t)¥ + F(t), 


A linear system is a system of differential equa- 


A non-homogeneous system of 
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ay fi Gil ++: Qin 
x= : ; f= ; ; A= 


In Fn QGm1 *** GAmn 
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Brine Tank Cascade 
Let brine tanks A, B, C' be given of volumes 20, 40, 60, respectively, as 


in Figure 1. 
water 


| 
[A =) 
B n>: 
—| 
v | Figure 1. Three brine 


tanks in cascade. 
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It is supposed that fluid enters tank A at rate r, drains from A to B 
at rate r, drains from B to C at rate r, then drains from tank C' at 
rate r. Hence the volumes of the tanks remain constant. Let r = 10, to 
illustrate the ideas. 


Uniform stirring of each tank is assumed, which implies uniform salt 
concentration throughout each tank. 


Let 21(t), v2(t), 73(t) denote the amount of salt at time ¢ in each tank. 
We suppose water containing no salt is added to tank A. Therefore, 
the salt in all the tanks is eventually lost from the drains. The cascade 
is modeled by the chemical balance law 


rate of change = inputrate — output rate. 


Application of the balance law, justified below in compartment analysis, 
results in the triangular differential system 


/ 1 

Ly = 75" 

ct ites, <a 
v9 = gut = qo? 
ee, el 
v3 = qv = 6°: 


The solution, to be justified later in this chapter, is given by the equations 


x(t) = a1 (O)e”, 
v(t) = —221(0)e~*/? + (x9(0) + 221(0))e~4, 


BOS : 24 (O)e~*/2 — 3(ar9(0) + 22 (0))e~*/4 


PGR Ope *71(0) + 3(wo(0) + 2a (0)))en*/®, 


Cascades and Compartment Analysis 


A linear cascade is a diagram of compartments in which input and 
output rates have been assigned from one or more different compart- 
ments. The diagram is a succinct way to summarize and document the 
various rates. 


The method of compartment analysis translates the diagram into a 
system of linear differential equations. The method has been used to 
derive applied models in diverse topics like ecology, chemistry, heating 
and cooling, kinetics, mechanics and electricity. 


The method. Refer to Figure 2. A compartment diagram consists of 
the following components. 


Variable Names Each compartment is labelled with a variable _X. 
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Arrows Each arrow is labelled with a flow rate R. 
Input Rate An arrowhead pointing at compartment X docu- 
ments input rate R. 


Output Rate An arrowhead pointing away from compartment X 
documents output rate R. 


0 2/2 
> 


Figure 2. Compartment 
analysis diagram. 


x 
: The diagram represents the 
13/6 classical brine tank problem of 
Figure 1. 


Assembly of the single linear differential equation for a diagram com- 
partment X is done by writing dX/dt for the left side of the differential 
equation and then algebraically adding the input and output rates to ob- 
tain the right side of the differential equation, according to the balance 
law 


dX : 
aE = sum of input rates — sum of output rates 


By convention, a compartment with no arriving arrowhead has input 
zero, and a compartment with no exiting arrowhead has output zero. 
Applying the balance law to Figure 2 gives one differential equation for 
each of the three compartments | 7; |, | v2 |, | x3 |. 


t= a 
1 2 1; 
oi. a 

Lg = 571 = 4°? 
Pe eee 

3 >= Av? = eo 


Recycled Brine Tank Cascade 


Let brine tanks A, B, C' be given of volumes 60, 30, 60, respectively, as 
in Figure 3. 


A Ey C 
| B Figure 3. Three brine tanks 
; -———_ in cascade with recycling. 


Suppose that fluid drains from tank A to B at rate r, drains from tank 
B to C at rate r, then drains from tank C to A at rate r. The tank 
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volumes remain constant due to constant recycling of fluid. For purposes 
of illustration, let r = 10. 


Uniform stirring of each tank is assumed, which implies uniform salt 


concentration throughout each tank. 


Let 21(t), v2(t), x3(t) denote the amount of salt at time ¢ in each tank. 
No salt is lost from the system, due to recycling. Using compartment 
analysis, the recycled cascade is modeled by the non-triangular system 


1 
a = er + Gs, 
‘ 1 1 
oe 6°! = 32 
1 
xv = ~X2 oS aes 
3 3 6 


The solution is given by the equations 

a(t) = cy + (cg — 2c3)e*/ cos(t/6) + (2c. + c3)e~/* sin(t/6), 
cy + (—2c2 — c3)e~*/3 cos(t/6) + (cz — 2c3)e~*/3 sin(t/6), 
13(t) = c1 + (co + 3c3)e~“/? cos(t/6) + (—3c2 + c3)e~“/? sin(t/6). 


At infinity, 71] = 73 = cj, 2 = c,/2. The meaning is that the total 
amount of salt is uniformly distributed in the tanks, in the ratio 2:1: 2. 


Pond Pollution 


Consider three ponds connected by streams, as in Figure 4. The first 
pond has a pollution source, which spreads via the connecting streams 
to the other ponds. The plan is to determine the amount of pollutant in 
each pond. 


Figure 4. Three ponds 1, 2, 3 
of volumes V,, V2, V3 connected 

f(t by streams. The pollution 
source f(t) is in pond 1. 


Assume the following. 


e Symbol f(t) is the pollutant flow rate into pond 1 (lb/min). 


e Symbols fi, f2, f3 denote the pollutant flow rates out of ponds 1, 
2, 3, respectively (gal/min). It is assumed that the pollutant is 
well-mixed in each pond. 
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e The three ponds have volumes Vj, V2, V3 (gal), which remain con- 
stant. 


e Symbols x1(t), v(t), x3(t) denote the amount (Ibs) of pollutant in 
ponds 1, 2, 3, respectively. 


The pollutant flux is the flow rate times the pollutant concentration, e.g., 
pond 1 is emptied with flux f; times 71(t)/Vi. A compartment analysis 
is summarized in the following diagram. 


F(t) fiti/Vi 
a LY SS x2 
Figure 5. Pond diagram. 
f3x3/V3 The compartment diagram 
/ fa2/Va represents the three-pond 
ws pollution problem of Figure 4. 


The diagram plus compartment analysis gives the following differential 
equations. 


w(t) = H a(t) Bn) +I) 
4) = Bn 2a, 
ie 2 a(t) — Boalt) 


For a specific numerical example, take f;/V; = 0.001, 1 < i < 3, and 
let f(t) = 0.125 lb/min for the first 48 hours (2880 minutes), thereafter 
f(t) =0. We expect due to uniform mixing that after a long time there 
will be (0.125)(2880) = 360 pounds of pollutant uniformly deposited, 
which is 120 pounds per pond. 


Initially, 71(0) = x2(0) = x3(0) = 0, if the ponds were pristine. The 
specialized problem for the first 48 hours is 


‘(t) = 0.001 x3(t) — 0.001 21 (t) + 0.125, 
x(t) = 0.001 2;(t) — 0.001 x(t), 

3(t) = 0.001 x(t) — 0.001 x3(t), 

1(0) x2(0) x3(0) 0. 


The solution to this system is 


125V3 3t 125 3t 125 t 
a(t) =a ( Fan (4 cos (45) ) EY 


9 2000 3 2000 3. 24? 
(i 250/3 ae V/3t Z t 
ne 9. \ 9000) 34? 


r3(t) = e~ 2000 & cos (5) ee an (He) + ae = Bee 


2000 9 2000 o4° 3B 


11.1 Examples of Systems 745 


After 48 hours elapse, the approximate pollutant amounts in pounds are 
x1(2880) = 162.30, x2(2880) = 119.61, 23(2880) = 78.08. 


It should be remarked that the system above is altered by replacing 0.125 
by zero, in order to predict the state of the ponds after 48 hours. The 
corresponding homogeneous system has an equilibrium solution 2(t) = 
xo(t) = x3(t) = 120. This constant solution is the limit at infinity of 
the solution to the homogeneous system, using the initial values 71 (0) = 
162.30, x2(0) © 119.61, x3(0) & 78.08. 


Home Heating 
Consider a typical home with attic, basement and insulated main floor. 


Attic 


Main 

Floor 
Figure 6. Typical home 

Basement with attic and basement. 
The below-grade basement 
and the attic are un-insulated. 
Only the main living area is 
insulated. 


It is usual to surround the main living area with insulation, but the attic 
area has walls and ceiling without insulation. The walls and floor in the 
basement are insulated by earth. The basement ceiling is insulated by 
air space in the joists, a layer of flooring on the main floor and a layer 
of drywall in the basement. We will analyze the changing temperatures 
in the three levels using Newton’s cooling law and the variables 


z(t) = Temperature in the attic, 
y(t) = Temperature in the main living area, 
x(t) = Temperature in the basement, 


t = Time in hours. 


Initial data. Assume it is winter time and the outside temperature 
in constantly 35°F during the day. Also assumed is a basement earth 
temperature of 45°F. Initially, the heat is off for several days. The initial 
values at noon (t = 0) are then (0) = 45, y(0) = z(0) = 35. 

Portable heater. A small electric heater is turned on at noon, with 
thermostat set for 100°F. When the heater is running, it provides a 20°F 
rise per hour, therefore it takes some time to reach 100°F (probably 
never!). Newton’s cooling law 
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Temperature rate = k(Temperature difference) 


will be applied to five boundary surfaces: (0) the basement walls and 
floor, (1) the basement ceiling, (2) the main floor walls, (3) the main 
floor ceiling, and (4) the attic walls and ceiling. Newton’s cooling law 
gives positive cooling constants ko, ki, ko, k3, k4 and the equations 


ve = ko(45-2) +k (y-2), 
yo = ky (a — y) + ko(35 — y) + ka(z — y) + 20, 
: a k3(y — z) + ka(35 — 2). 


The insulation constants will be defined as kg = 1/2, ky = 1/2, ko = 1/4, 
k3 = 1/4, kg = 3/4 to reflect insulation quality. The reciprocal 1/k 
is approximately the amount of time in hours required for 63% of the 
temperature difference to be exchanged. For instance, 4 hours elapse for 
the main floor. The model: 


1 1 

go = ae oe £), 

yo = 5 (t —Y) + 7 (35 — y) + F(z —y) + 20, 
1 

a= racy z)4 (85 2) 


The homogeneous solution in vector form is given in terms of constants 
a=-1+4+¥75/4,b=1- 5/4, and arbitrary constants c,, c2, c3 by the 
formula 


rn(t) —1 2 2 
yp(t) | =cae* 0 | +eoe% | WS | tee | —/5 
zp(t) 2 1 1 


p(t) geo 
Yp(t) = ne 
Zp(t) Se 


The homogeneous solution has limit zero at infinity, hence the temper- 
atures of the three spaces hover around x = 56.4, y = 67.7, z = 43.2 
degrees Fahrenheit. Specific information can be gathered by solving for 
C1, C2, cz according to the initial data x(0) = 45, y(0) = z(0) = 35. The 
answers are 


D5. 2 
Cl , © 7 + ¥5: C3 5 


Underpowered heater. To the main floor each hour is added 20°F, but 
the heat escapes at a substantial rate, so that after one hour y = 68°F. 
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After five hours, y © 68°F. The heater in this example is so inadequate 
that even after many hours, the main living area is still under 69°F. 


Forced air furnace. Replacing the space heater by a normal furnace 
adds the difficulty of switches in the input, namely, the thermostat 
turns off the furnace when the main floor temperature reaches 70°F, 
and it turns it on again after a 4°F temperature drop. We will suppose 
that the furnace has four times the BTU rating of the space heater, 
which translates to an 80°F temperature rise per hour. The study of 
the forced air furnace requires two differential equations, one with 20 
replaced by 80 (DE 1, furnace on) and the other with 20 replaced by 0 
(DE 2, furnace off). The plan is to use the first differential equation on 
time interval 0 < t < t,, then switch to the second differential equation 
for time interval t] < t < tg. The time intervals are selected so that 
y(t) = 70 (the thermostat setting) and y(t2) = 66 (thermostat setting 
less 4 degrees). Numerical work gives the following results. 


Time in minutes | Main floor temperature | Model | Furnace 


31.6 70 DE 1 on 
40.9 66 DE 2 off 
45.3 70 DE 1 on 
54.6 66 DE 2 off 


The reason for the non-uniform times between furnace cycles can be 
seen from the model. Each time the furnace cycles, heat enters the main 
floor, then escapes through the other two levels. Consequently, the initial 
conditions on each floor applied to models 1 and 2 are changing, resulting 
in different solutions to the models on each switch. 


Chemostats and Microorganism Culturing 


A vessel into which nutrients are pumped, to feed a microorganism, 
is called a chemostat!. Uniform distributions of microorganisms and 
nutrients are assumed, for example, due to stirring effects. The pumping 
is matched by draining to keep the volume constant. 


'The October 14, 2004 issue of the journal Nature featured a study of the co- 
evolution of a common type of bacteria, Escherichia coli, and a virus that infects 
it, called bacteriophage T7. Postdoctoral researcher Samantha Forde set up ” micro- 
bial communities of bacteria and viruses with different nutrient levels in a series of 
chemostats — glass culture tubes that provide nutrients and oxygen and siphon off 
wastes.” 
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Input Feed Output Effluent 


Figure 7. A Basic 
Chemostat. A stirred 
bio-reactor operated as a 
chemostat, with continuous inflow 
and outflow. The flow rates are 
controlled to maintain a constant 
culture volume. 


In a typical chemostat, one nutrient is kept in short supply while all 
others are abundant. We consider here the question of survival of the 
organism subject to the limited resource. The problem is quantified as 
follows: 


x(t) = the concentration of the limited nutrient in the vessel, 


y(t) = the concentration of organisms in the vessel. 


A special case of the derivation in J.M. Cushing’s text for the organism 
E. Coli? is the set of nonlinear differential equations® 


1 
x’ = —0.075x + (0.075)(0.005) — —g(zx)y, 


(2) ; 63 
y = —0.075y + g(x)y, 


where g(x) = 0.682(0.0016 + x)~!. Of special interest to the study of 
this equation are two linearized equations at equilibria, given by 


(3) ui, = —0.075 uz — 0.008177008175 ue, 
ul, = 0.4401515152 ua, 


—1.690372243 v; — 0.001190476190 vo, 
101.7684513 v1. 


vj 
Vg 


(4) 


Although we cannot solve the nonlinear system explicitly, nevertheless 
there are explicit formulas for u,, ug, v1, v2 that complete the picture of 
how solutions x(t), y(t) behave at t = oo. The result of the analysis is 
that E. Coli survives indefinitely in this vessel at concentration y ~ 0.3. 


?In a biology Master’s thesis, two strains of Escherichia coli were grown in a glucose- 
limited chemostat coupled to a modified Robbins device containing plugs of silicone 
rubber urinary catheter material. Reference: Jennifer L. Adams and Robert J. C. 
McLean, Applied and Environmental Microbiology, September 1999, p. 4285-4287, 
Vol. 65, No. 9. 

3More details can be found in The Theory of the Chemostat Dynamics of Microbial 
Competition, ISBN-13: 9780521067348, by Hal Smith and Paul Waltman, June 2008. 


11.1 Examples of Systems 749 


air inlet 


Feed Reservoir 


pump 
air inlet 
< t~ heater /cooler 


stirring bar 


Culture vessel 


overflow 


magnetic stirrer 


Effluent reservoir 


Figure 8. Laboratory Chemostat. 
The components are the Feed reservoir, which contains the nutrients, a stirred 


chemical reactor labeled the Culture vessel, and the Effluent reservoir, 
which holds the effluent overflow from the reactor. 


Irregular Heartbeats and Lidocaine 


The human malady of ventricular arrhythmia or irregular heartbeat 
is treated clinically using the drug lidocaine. 


" Xylocaine 


0.5% {lidocaine 
0/0 Injocti SF 


Figure 9. Xylocaine label, a brand name for 
the drug lidocaine. 


To be effective, the drug has to be maintained at a bloodstream concen- 
tration of 1.5 milligrams per liter, but concentrations above 6 milligrams 
per liter are considered lethal in some patients. The actual dosage de- 
pends upon body weight. The adult dosage maximum for ventricular 
tachycardia is reported at 3 mg/kg.* The drug is supplied in 0.5%, 1% 
and 2% solutions, which are stored at room temperature. 


A differential equation model for the dynamics of the drug therapy uses 


x(t) = amount of lidocaine in the bloodstream, 


y(t) = amount of lidocaine in body tissue. 


A typical set of equations, valid for a special body weight only, appears 
below; for more detail see J.M. Cushing’s text [Cushing (2004)]. 
(5) x(t) = —0.09x(t) + 0.038y(t), 

y'(t) = 0.066x(t) — 0.038y(t). 


‘Source: Family Practice Notebook, http://www.fpnotebook.com/. The au- 
thor is Scott Moses, MD, who practises in Lino Lakes, Minnesota. 
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The physically significant initial data is zero drug in the bloodstream 
x(0) = 0 and injection dosage y(0) = yo. The answers: 


a(t) = —0.3367ype~ 9-170 +. 0.3367 ype7 0-007, 
y(t) = 0.2696ype~ 0-124 + 0.7304 ype 0007". 


The answers can be used to estimate the maximum possible safe dosage 
yo and the duration of time that the drug lidocaine is effective. 


Nutrient Flow in an Aquarium 


Consider a vessel of water containing a radioactive isotope, to be used as 
a tracer for the food chain, which consists of aquatic plankton varieties 
A and B. 


Plankton are aquatic organisms that drift with the currents, typically 
in an environment like Chesapeake Bay. Plankton can be divided into 
two groups, phytoplankton and zooplankton. The phytoplankton are 
plant-like drifters: diatoms and other alga. Zooplankton are animal-like 
drifters: copepods, larvae, and small crustaceans. 


Figure 10. Left: Bacillaria 
paxillifera, phytoplankton. 
Right: Anomura Galathea 
zoea, zooplankton. 


Let 


x(t) = isotope concentration in the water, 
y(t) = isotope concentration in A, 


z(t) = isotope concentration in B. 


Typical differential equations are 


x’ (t) = —32(t) + 6y(t) + 52(0), 
y(t) = 2a(t) — 12y(t), 
z'(t) = x(t) + 6y(t) — 5z(t). 


The answers are 


x(t) = 6c, + (1+ V6)ege1+ V9)" + (1 — VB)ege(—1O-VOt, 
y(t) =cr+ cge(— 10+ v6) + ege(—10-V6)t 


12 
2(t) = ae = (2 +v 15) cge(—10+ v6) + (-2 + V15) cge(—10-V6)t_ 
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The constants cj, c2, cg are related to the initial radioactive isotope 
concentrations x(0) = 20, y(0) = 0, z(0) = 0, by the 3 x 3 system of 
linear algebraic equations 


6c, + (1+ V6)co + (1—V6)ce3 = 20, 
Gi. = cg + cz = 0, 
12 


=a = (2+V15)q + (-2+V15)cs 


I 
S 


Biomass Transfer 


Consider a European forest having one or two varieties of trees. We 
select some of the oldest trees, those expected to die off in the next few 
years, then follow the cycle of living trees into dead trees. The dead trees 
eventually decay and fall from seasonal and biological events. Finally, 
the fallen trees become humus. 


Figure 11. Forest Biomass. Total biomass is a parameter used to assess 
atmospheric carbon that is harvested by trees. Forest management uses biomass 
subclasses to classify fire risk. 


Let variables x, y, z, t be defined by 


x(t) = biomass decayed into humus, 
y(t) = biomass of dead trees, 
z(t) = biomass of living trees, 


t = time in decades (decade = 10 years). 


A typical biological model is 


x(t) = —a(t) + 3y(t), 
y'(t) = —8y(t) + 52(t), 
z2'(t) = —52(t). 
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Suppose there are no dead trees and no humus at t = 0, with initially zo 
units of living tree biomass. These assumptions imply initial conditions 
x(0) = y(0) = 0, z(0) = z. The solution is 


15 
i) = 3 (c of eae " 
5 
y(t) = 52 (-e* +e), 
z(t) = ae". 


The live tree biomass z(t) = ze >! decreases according to a Malthusian 
decay law from its initial size zg. It decays to 60% of its original biomass 
in one year. Interesting calculations that can be made from the other 
formulas include the future dates when the dead tree biomass and the 
humus biomass are maximum. The predicted dates are approximately 
2.5 and 8 years hence, respectively. 


The predictions made by this model are trends extrapolated from rate 
observations in the forest. Like weather prediction, it is a calculated 
guess that disappoints on a given day and from the outset has no pre- 
dictable answer. 


Total biomass is considered an important parameter to assess atmo- 
spheric carbon that is harvested by trees. Biomass estimates for forests 
since 1980 have been made by satellite remote sensing data with instances 
of 90% accuracy (Science 87(5), September 2004). 

Pesticides in Soil and Trees 


A Washington cherry orchard was sprayed with pesticides. 


Figure 12. Cherries in June. 


Assume that a negligible amount of pesticide was sprayed on the soil. 
Pesticide applied to the trees has a certain outflow rate to the soil, and 
conversely, pesticide in the soil has a certain uptake rate into the trees. 
Repeated applications of the pesticide are required to control the insects, 
which implies the pesticide levels in the trees varies with time. Quantize 
the pesticide spraying as follows. 


11.1 Examples of Systems 753 


x(t) = amount of pesticide in the trees, 
y(t) = amount of pesticide in the soil, 
r(t) = amount of pesticide applied to the trees, 


¢t = time in years. 


A typical model is obtained from input-output analysis, similar to the 
brine tank models: 


(t) 


vid 2x(t) — y(t) + r(t), 
y'(t) 


= 2x(t) — 3y(t). 

In a pristine orchard, the initial data is 7(0) = 0, y(0) = 0, because the 
trees and the soil initially harbor no pesticide. The solution of the model 
obviously depends on r(t). The nonhomogeneous dependence is treated 
by the method of variation of parameters infra. Approximate formulas 
are 


t 

a(t) | (1.10¢160—) - 0.12¢-?-4-») | r(u)du, 
0 
t 

y(t) = | (0.49e1-60-™) — 0.49¢e~?-6(t-») ) r(u)du. 
0 


The exponential rates 1.6 and —2.6 represent respectively the accumu- 
lation of the pesticide into the soil and the decay of the pesticide from 
the trees. The application rate r(t) is typically a step function equal to 
a positive constant over a small interval of time and zero elsewhere, or a 
sum of such functions, representing periodic applications of pesticide. 


Forecasting Prices 


A cosmetics manufacturer has a marketing policy based upon the price 
x(t) of its salon shampoo. 


Figure 13. Salon 
shampoo sample. 

The marketing strategy for 
the shampoo is to set the 
price z(t) dynamically to 
reflect demand for the 
product. 
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The production P(t) and the sales S(t) are given in terms of the price 
x(t) and the change in price z'(t) by the equations 


P(t) =4- “a(t) — 82'(t) (Production), 


S(t) = 15 — 4x(t) — 22'(t) (Sales). 
The differential equations for the price x(t) and inventory level I(t) are 
x(t) = k(I(t) — Ip), 
I'(t) = P(t) — S(t). 


The inventory level J) = 50 represents the desired level. The equations 
can be written in terms of x(t), I(t) as follows. 


et) = kI(t) — klp, 
roy = x(t) GRR) ae Bhp Al. 


If k = 1, «(0) = 10 and J(0) = 7, then the solution is given by 


= 44 86 —13t/2 
WO Hag tage 
I(t) = 50 — 43¢783#/2, 


The forecast of price x(t) + 3.39 dollars at inventory level I(t) + 50 is 
based upon the two limits 


too t-0o 


ioe - lim I(t) = 50. 


Coupled Spring-Mass Systems 


Three masses are attached to each other by four springs as in Figure 14. 
Figure 14. Three masses 


ky kag kg ig : 
Ts A SUS | SON connected by springs. The masses 
slide along a frictionless horizontal 
mi, mg m3 surface. 


The analysis uses the following constants, variables and assumptions. 


Mass The masses m1, M2, m3 are assumed to be point masses 
Constants concentrated at their center of gravity. 

Spring The mass of each spring is negligible. The springs op- 
Constants erate according to Hooke’s law: Force = k(elongation). 


Constants ki, ko, k3, k4 denote the Hooke's constants. 

The springs restore after compression and extension. 
Position The symbols x(t), x2(t), v3(t) denote the mass posi- 
Variables tions along the horizontal surface, measured from their 

equilibrium positions, plus right and minus left. 
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Fixed Ends _ The first and last spring are attached to fixed walls. 


The competition method is used to derive the equations of motion. 
In this case, the law is 


Newton’s Second Law Force = Sum of the Hooke’s Forces. 


The model equations are 


may (t) = —ky21(t) + ke[wo(t) — 21 (0), 
(6) mgxg(t) = —ke[xo(t) — x1(t)] + ka[xs(t) — x2(t)], 
m3x5(t) = —keg [x3(t) = x(t)] = kga3(t). 


The equations are justified in the case of all positive variables by observ- 
ing that the first three springs are elongated by 21, r2 — 11, ©3 — 22, 
respectively. The last spring is compressed by x3, which accounts for the 
minus sign. 


Another way to justify the equations is through mirror-image symmetry: 
interchange ky <—> ky, ko <> kg, x1 <> 23, then equation 2 should be 
unchanged and equation 3 should become equation 1. 


Matrix Formulation. System (6) can be written as a second order 
vector-matrix system 


My] 0 O xt —ky = ko ko 0 Ly 
0 mg 0 oe = ko —ko = kg k3 Hop) 
0 0 M3, ae 0 kg —kg = ka x3 


More succinctly, the system is written as 
Mx" (t) = KX (t) 


where the displacement X, mass matrix M and stiffness matrix K 
are defined by the formulas 


Ly my 0 0 —ky —_ ko ko 0 
X= TQ), M= 0 mg 0 5 K= ko —ky — kg k3 
£3 0 0 ms 0 kg = —k3 — kg 


Numerical example. Let m, = 1, m2 = 1, m3 = 1, ky = 2, ko = 1, 
k3 = 1, kg = 2. Then the system is given by 


et -3 1 O Ly 
ot |={ 1-2 1 |] 2 
xy 0 1 -8 r3 
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The vector solution is given by the formula 


LY 1 
x2 = (a,cost+6,sint) | 2 
£3 1 


1 
+ (a2 cos V3t + b2 sin v3t) 0 
—l 


+ (a3 cos 2¢ + bgsin2t) | —1 |, 
where aj, a2, a3, 61, bo, b3 are arbitrary constants. 


Boxcars 


A special case of the coupled spring-mass system is three boxcars on a 
level track connected by springs, as in Figure 15. 


k k Figure 15. Three identical 
aaa, i si Sa boxcars connected by 
m m m identical springs. 


Except for the springs on fixed ends, this problem is the same as the one 
of the preceding illustration, therefore taking kj = ky = 0, ko = k3 =k, 
my, = M2 = m3 = mM gives the system 


m0 0 xy —k k 0 Ty 
OmO|{as}=] k-2k k| | 2 
00m} \a3 0 k-k} \a3 
Take k/m = 1 to obtain the illustration 
—1 1 0 
— 1 -2 1] x, 
0 1 -1 
which has vector solution 
1 1 
X = (a,+6)t)| 1 | + (agcost + besint) 0 
1 —1 
1 
+ (a3 cos V/3t + b3 sin v3t) —2 |, 
1 


where ay, a2, a3, 61, bo, b3 are arbitrary constants. 
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The solution expression can be used to discover what happens to the 
boxcars when the springs act normally upon compression but disengage 
upon expansion. An interesting physical situation is when one car moves 
along the track, contacts two stationary cars, then transfers its momen- 
tum to the other cars, followed by disengagement. 


Monatomic Crystals 


Figure 16. A Crystal Model. 


The n crystals are identical masses m assumed connected by equal springs of 
Hooke’s constant k. The last mass is connected to the first mass. 


The scalar differential equations for Figure 16 are written for mass po- 
sitions 11,...,0%p, with 79 = In, Xn41 = X1 to account for the ring of 
identical masses (periodic boundary condition). Then for k = 1,...,n 
dx}, 
are 
These equations represent a system x” = Ax, where the symmetric ma- 
trix of coefficients A = M~!K is given for n = 5 and k/m =1 by 


—2 1 O 0 


= k(te41 — Be) + h(ee_-1 — By) = (eR-1 — 20% + Fe41). 


1 

1 —-2 1 0 0 

A= 0 1 -—2 1 0 

0 0 1 —2 1 

1 0 0 1 —2 
—2 1 i 

Ifn = 3 and k/m = 1, then A= 1—2 1] and the solutions x1, x2, 

1 1-2 


x3 are linear combinations of the functions 1, t, cos V3t, sin V3¢. 


Electrical Network I 


Consider the L.R-network of Figure 17. 


Figure 17. An 
electrical network. 
There are three 


a4 Ty 


Ry 


resistors R,, Ro, Rg 
and three inductors 
Ih, La, D3. The 
currents 71, 22, 723 are 
defined between 
nodes (black dots). 
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The derivation of the differential equations for the loop currents 71, 72, 
i3 uses Kirchhoff’s laws and the voltage drop formulas for resistors and 
inductors. The black dots in the diagram are the nodes that determine 
the beginning and end of each of the currents 71, 12, 13. Currents are 
defined only on the outer boundary of the network. Kirchhoff’s node law 
determines the currents across L2, L3 (arrowhead right) as ig — 71 and 
iz —i1, respectively. Similarly, i2—73 is the current across R; (arrowhead 
down). Using Ohm’s law Vg = RI and Faraday’s law Vp = LI’ plus 
Kirchhoff’s loop law algebraic sum of the voltage drops is zero around a 
closed loop (see the maple code below), we arrive at the model 


ne (BeBae (BB) 
3 Ey. fee * te Tage 
Wa (> =)i (> +p) 
: Lz iy) * is Dj ig)" 
A computer algebra system is helpful to obtain the differential equations 


from the closed loop formulas. Part of the theory is that the number of 
equations equals the number of holes in the network, called the connec- 


tivity. Here’s some maple code for determining the equations in scalar 
and also in vector-matrix form. 


loop1:=L1*D(i1)+R3*i3+R2*i2=0; 

loop2:=L2*D (i2)-L2*D(i1)+R1* (i2-i3) +R2*i2=0; 
lLoop3:=L3*D (i3) -L3*D (11) +R3*i3+R1* (i3-i2)=0; 
f1:=solve(loop1,D(i1)); 
£2:=solve(subs(D(i1)=f1,loop2) ,D(i2)); 
£3:=solve(subs (D(i1)=f1,loop3) ,D(i3)); 
with(linalg): 
jacobian([f1,f£2,£3] , [i1,i2,i3]); 


Electrical Network II 


Consider the L.R-network of Figure 18. This network produces only two 
differential equations, even though there are three holes (connectivity 
3). The derivation of the differential equations parallels the previous 
network, so nothing will be repeated here. 


A computer algebra system is used to obtain the differential equations 
from the closed loop formulas. Below is maple code to generate the 
equations 7, = fi, 14 = fo, 13 = fs. 


loop1:=L1*D(i1)+R2* (i1-i2)+R1*(i1-i3)=0; 
loop2:=L2*D (i2)+R3* (i2-i3)+R2* (i2-i1)=0; 
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loop3 : =R3* (i3-i2)+R1* (i3-i1)=E; 
£3:=solve(loop3,i3) ; 
f£1:=solve (subs (i3=£3,loop1) ,D(i1)); 
£2:=solve (subs (i3=£3,loop2) ,D(i2)); 


: OST 
—S ay Ln 
Ry Ro 
—_—$__ WV 
©) As " 
23<=———_ [ 12 = 


Figure 18. An electrical network. 
There are three resistors R,, Ro, R3, two inductors L,, Lz and a battery E. 
The currents i1, 72, 13 are defined between nodes (black dots). 


The model, in the special case DL; = Lz = 1 and Ri = Ro = R3 = R: 


y= = an, + oy + a 
| 2 1 2 2” 
re 
> = 2 1 9 2 2” 

= i. 4 ‘ 
eo fl 92 aR’ 


It is easily justified that the solution of the differential equations for 
initial conditions 71(0) = 72(0) = 0 is given by 


= Ft, ig(t) = Ze 


Logging Timber by Helicopter 


Certain sections of National Forest in the USA do not have logging ac- 
cess roads. In order to log the timber in these areas, helicopters are 
employed to move the felled trees to a nearby loading area, where they 
are transported by truck to the mill. The felled trees are slung beneath 
the helicopter on cables. 


Figure 19. Helicopter logging. 
Left: An Erickson helicopter lifts felled 
trees. 

Right: Two trees are attached to the 
cable to lower transportation costs. 
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The payload for two trees approximates a double pendulum, which oscil- 
lates during flight. The angles of oscillation 61, 62 of the two connecting 
cables, measured from the gravity vector direction, satisfy the following 
differential equations, in which g is the gravitation constant, m1, ma 
denote the masses of the two trees and Lj, Lz are the cable lengths. 


(m4 + m2) L204 iad mL L264 + (m1 + m2)L4g9, = 0, 
mL L26f =F m2L264 + m2L2g02 = (); 


This model is derived assuming small displacements 6), 62, that is, 
sin @ = @ for both angles, using the following diagram. 


01 
Ly Figure 20. Logging Timber by Helicopter. 

mo Bee The cables have lengths L,, Lz. The angles @;, 02 are 
: A» measured from vertical. 


The lengths £1, Lz are adjusted on each trip for the length of the trees, 
so that the trees do not collide in flight with each other nor with the 
helicopter. Sometimes, three or more smaller trees are bundled together 
in a package, which is treated here as identical to a single, very thick 
tree hanging on the cable. 

Vector-matrix model. The angles 4;, 92 satisfy the second-order 
vector-matrix equation 


(mi + m2)L1 meL2 0, \" _ _{ mg+meg 0 a 
Ty Le A a 0 g 02 


This system is equivalent to the second-order system 


7 mig + mag m2g 
0, = Lym4 Lim, 6; 
i | mig+meg (m+ m2) g tal 
Lym, Lym, 


Earthquake Effects on Buildings 


A horizontal earthquake oscillation F(t) = Fo cos wt affects each floor of 
a 5-floor building; see Figure 21. The effect of the earthquake depends 
upon the natural frequencies of oscillation of the floors. 


In the case of a single-floor building, the center-of-mass position x(t) 
of the building satisfies ma” + kx = E and the natural frequency of 
oscillation is \/k/m. The earthquake force FE is given by Newton’s second 
law: E(t) = —mF"(t). Ifw & \/k/m, then the amplitude of x(t) is large 
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compared to the amplitude of the force E. The amplitude increase in 
x(t) means that a small-amplitude earthquake wave can resonant with 
the building and possibly demolish the structure. 


Figure 21. A 5-Floor 
Building. 

A horizontal earthquake wave F 
affects every floor. The actual wave 
has wavelength many times larger 
=> F than the illustration. 


PN Ww -F ol 


[co | 


The following assumptions and symbols are used to quantize the oscilla- 
tion of the 5-floor building. 


e Each floor is considered a point mass located at its center-of-mass. 
The floors have masses ™m} , ..., ™5. 


e Each floor is restored to its equilibrium position by a linear restor- 
ing force or Hooke’s force —k(elongation). The Hooke’s constants 


are ky, tees ks. 
e The locations of masses representing the 5 floors are 11, ..., Z5. 
The equilibrium position is 7] =--- = 25 = 0. 


e Damping effects of the floors are ignored. This is a frictionless 
system. 


The differential equations for the model are obtained by competition: 
the Newton’s second law force is set equal to the sum of the Hooke’s 
forces and the external force due to the earthquake wave. This results in 
the following system, where kg = 0, Ej = mjF" for j = 1,2,3,4,5 and 
F = Focoswt. 


myx —(ky + k)a4 + kore + FE}, 


mot, = koxy — (ko +kz)xe + kgx3 + Eo, 
m3t3 = kgx — (kg + ka)a3 + kava t+ Es, 
maxg = kavg — (kat ks)aa + kexts + Ea, 


mst, = ksxq — (ks + ke)a5 + Es. 


In particular, the equations for a floor depend only upon the neighboring 
floors. The bottom floor and the top floor are exceptions: they have just 
one neighboring floor. 


Vector-matrix second order system. Define 


My 0 0 0 0 XY Ey 

0 mg 0 0 0 x2 E» 
M=]!] 0 0 mz, 0 O |, X=] 23 |, H=] B; |, 

0 0 0 M4 0 LA Ey 

0 0 0 0 m5, XS Es 
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—k, — ke ko 0 0 0 
ko ay ee ks 0 0 
oe 0 kg = hs = hy ka 0 
0 0 ka —ka — ks ks 
0 0 0 ks he ky 


In the last row, kg = 0, to reflect the absence of a floor above the fifth. 
The second order system is 


Mx" (t) = KX(t) + H(t). 


The matrix M is called the mass matrix and the matrix Kk is called 
the Hooke’s matrix. The external force H(t) can be written as a 
scalar function E(t) = —F’”(t) times a constant vector: 


H(t) = —w*Focoswt | m3 


Identical floors. Let us assume that all floors have the same mass 
m and the same Hooke’s constant k. Then M = mI and the equation 
becomes 


—2k; k 0 0 0 1 

k —2k k 0 60 it 

"=m 0 k —2k k 0 | ¥ — Fow? cos(wt) | 1 
0 0 k —2k k; 1 

0 0 0 k —-k 1 


The Hooke’s matrix K is symmetric (K’ = K) with negative entries 
only on the diagonal. The last diagonal entry is —k (a common error is 
to write —2k). 


Particular solution. The method of undetermined coefficients predicts 
a trial solution x,(t) = € coswt, because each differential equation has 
nonhomogeneous term —Fow* coswt. The constant vector € is found 
by trial solution substitution. Cancel the common factor coswt in the 
substituted equation to obtain the equation (m7! K 4+ w? I) € = Fowb, 
where b is column vector of ones in the preceding display. Let B (w) = 


dj(B 
mK +w?TI. Then the formula B~! = rin gives 
. adj(B(w)) - 
= Fow b. 
Oe det (B(w)) 


The constant vector € can have a large magnitude when det(B(w)) = 0. 
This occurs when —w? is nearly an eigenvalue of m7!K. 
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Homogeneous solution. The theory of this chapter gives the homo- 
geneous solution X;(t) as the sum 


5 
cht) = So (a; cos wt + b; sinw;t)V ; 
j=l 


where r = w; and V =v; 4 0 satisfy 


and the values w1, ..., ws are found by solving the determinant equation 
det((1/m)K + wI) = 0, to obtain the values in Table 1. 


Table 1. The natural frequencies for the special case k/m = 10. 


Frequency Value 
Ww 0.900078068 
We 2.627315231 
W3 4.141702938 
Ww 5.320554507 
Ws 6.068366391 


General solution. Superposition implies x(t) = X,(t) + X,(t). Both 
terms of the general solution represent bounded oscillations. 


Resonance effects. The special solution X,(t) can be used to ob- 
tain some insight into practical resonance effects between the incoming 
earthquake wave and the building floors. When w is close to one of the 
frequencies w,, ..., W5, then the amplitude of a component of xX, can 
be very large, causing the floor to take an excursion that is too large to 
maintain the structural integrity of the floor. 


The physical interpretation is that an earthquake wave of the proper 
frequency, having time duration sufficiently long, can demolish a floor 
and hence demolish the entire building. The amplitude of the earthquake 
wave does not have to be large: a fraction of a centimeter might be 
enough to start the oscillation of the floors. 
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Earthquakes and Tsunamis 


Seismic wave shape was studied for first order equations in Chapter 2, 
page 151. Recorded here are some historical notes about seismic waves 
and earthquake events. 


The original Richter scale, with deprecated use in seismology, was 
invented by seismologist C. Richter to rank earthquake power. 


The moment magnitude scale (My) has largely replaced the original 
Richter scale and its modified versions. The highest reported magnitude 
is 9.5 Mw by the United States Geological Survey for the Concepci6én, 
Chile earthquake of May 22, 1960. News reports and the general public 
still refer to earthquake magnitude using the term Richter Scale. 


The Sumatra earthquake of December 26, 2004 occurred close to a deep- 
sea trench, a subduction zone where one tectonic plate slips beneath 
another. Most of the earthquake energy is released in these areas as the 
two plates grind towards each other. Estimates of magnitude 8.8 My 
to 9.3 My followed the event. The US Geological Survey estimated 
9.2 My. 


The Chile earthquake and tsunami of 1960 has been documented well. 
Here is an account by Dr. Gerard Fryer of the Hawaii Institute of Geo- 
physics and Planetology, in Honolulu. 


The tsunami was generated by the Chile earthquake of May 22, 
1960, the largest earthquake ever recorded: it was magnitude 9.6. 
What happened in the earthquake was that a piece of the Pacific 
seafloor (or strictly speaking, the Nazca Plate) about the size of 
California slid fifty feet beneath the continent of South America. 
Like a spring, the lower slopes of the South American continent 
offshore snapped upwards as much as twenty feet while land along 
the Chile coast dropped about ten feet. This change in the shape of 
the ocean bottom changed the shape of the sea surface. Since the 
sea surface likes to be flat, the pile of excess water at the surface 
collapsed to create a series of waves — the tsunami. 


The tsunami, together with the coastal subsidence and flooding, 
caused tremendous damage along the Chile coast, where about 
2,000 people died. The waves spread outwards across the Pa- 
cific. About 15 hours later the waves flooded Hilo, on the island 
of Hawaii, where they built up to 30 feet and caused 61 deaths 
along the waterfront. Seven hours after that, 22 hours after the 
earthquake, the waves flooded the coastline of Japan where 10-foot 
waves caused 200 deaths. The waves also caused damage in the 
Marquesas, in Samoa, and in New Zealand. Tidal gauges through- 
out the Pacific measured anomalous oscillations for about three 
days as the waves bounced from one side of the ocean to the other. 
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11.2 Basic First-order System Methods 


Solving 2 x 2 Systems 


It is shown here that any constant linear system 


ae _f ab 
u = Au, A=(2 5) 


can be solved by one of the following elementary methods. 


(a) The integrating factor method for y/ = p(x)y + q(z). 


(b) The second order constant coefficient formulas in Theo- 
rem 45, Chapter 5. 


Triangular A. Let’s assume b = 0, so that A is lower triangular. The 
upper triangular case is handled similarly. Then ti’ = Ati has the scalar 
form 

uy, = auy, 

us = cuit dug. 


The first differential equation is solved by the growth/decay formula: 


ui(t) = uge™. 


Then substitute the answer just found into the second differential equa- 
tion to give 

ub = dug + cuge™. 
This is a linear first order equation of the form y' = p(x)y + q(x), to be 
solved by the integrating factor method. Therefore, a triangular system 
can always be solved by the first order integrating factor method. 


An illustration. Let us solve i’ = AU for the triangular matrix 


1 0 ; Uy = th, 
A=( 3 a representing ee SS aoe 


The first equation u = uz has solution uy, = cje’. The second equation 
ul, = 2u1 + uz becomes upon substitution of uj = ce’ the new equation 


/ 
Uy = Ieie" + ua, 


which is a first order linear differential equation with linear integrating 
factor method solution ug = (2cit + ca)et. The general solution of ti’ = 
Au in scalar form is 


uit = cet, uo = 2c\ te’ + cet. 
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The vector form of the general solution is 


aa=a( su )re( 2). 


The vector basis is the set 


s{(s0).(3)} 


Non-Triangular A. In order that A be non-triangular, both b 4 0 
and c # 0 must be satisfied. The scalar form of the system a’ = Ad is 


u, = au,t bua, 
uy = curt dug. 


Theorem 1 (Solving Non-Triangular ui’ = At) 
Solutions ui, ug of i’ = Aw are linear combinations of the list of Euler 
solution atoms obtained from the roots r of the quadratic equation 


det(A — rI) = 0. 


Proof: The method: differentiate the first equation, then use the equations to 
eliminate uz, u5. The result is a second order differential equation for u;. The 
same differential equation is satisfied also for wz. The details: 


ul = au), + bus Differentiate the first equation. 
= au, + beu, + bdug Use equation uy = cu, + dug. 
= au, + beu, + d(u, — au) Use equation ui = au, + bua. 
=(a+d)u + (bc— ad)uy Second order equation for u; found 


The characteristic equation of ui — (a + d)u}, + (ad — bc)u1 = 0 is 
r? —(a+d)r + (be— ad) =0. 


Finally, we show the expansion of det(A — rJ) is the same characteristic poly- 
nomial: 
a-rT b 
det(A—rI) = tae 
= (a—r)(d—r) — bce 
= r*—(at+d)r+ad-— be. 


The proof is complete. 
The reader can verify that the differential equation for uz or ug is exactly 
u” — trace(A)u’ + det(A)u = 0. 


Assume below that A is non-triangular, meaning b #4 0 and c £ 0. 


Finding u,. Apply the second order formulas, Theorem 45 in Chapter 
5, to solve for uz. This involves writing a list of Euler solution atoms 
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corresponding to the two roots of the characteristic equation r? — (a + 
d)r + ad — bc = 0, followed by expressing ui as a linear combination of 
the two Euler atoms. 


Finding wz. Isolate u2 in the first differential equation by division: 


ug = —(u} — aur). 


ole 


The two formulas for u1, ug represent the general solution of the system 
u’ = Au, when A is 2 x 2. 


An illustration. Let’s solve i’ = AW when 


A= ae representin ty = wa + 2a, 
~~) 2 1 iy? ©P 8 us = 2u,+ us. 


The equation det(A — rl) = 0 is (1—r)? —4 =0 with roots r = —1 and 
r = 3. The Euler solution atom list is L = {e~',e*"}. Then the linear 
combination of Euler atoms is uy = cye~’ + ce". The first equation 
ui = uy; + 2uz2 implies ug = 4(u, — ui). The scalar general solution of 
u’ = Adi is then 


Uy = ce ¢ + coe", uz = —cje* + éer". 


In vector form, the general solution is 


“ et et 
u=C, et + C2 est . 


Triangular Methods 


Diagonal nxn matrix A = diag(a),...,a,,). Then the system X/ = Ax 
is a set of uncoupled scalar growth/decay equations: 


v(t) = aixi(t), 
v(t) = azxa(t), 
ee) = AnUn(t). 


The solution to the system is given by the formulas 


a(t) = cye™", 
xo(t) = cge%', 
Balt) = eget". 


The numbers cj, ..., Cn, are arbitrary constants. 
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Triangular n x n matrix A. Ifa linear system ¥’ = AX has a square 
triangular matrix A, then the system can be solved by first order scalar 
methods. To illustrate the ideas, consider the 3 x 3 linear system 


2 
pian Bs em em a ea 

444 
The coefficient matrix A is lower triangular. In scalar form, the system 
is given by the equations 


w(t) = 2a (t), 


v(t) = 321(t) + 3za(t), 
£3(t) = 421 (t) + 4x9 (t) + 4x3(t). 


A recursive method. The system is solved recursively by first order 
scalar methods only, starting with the first equation x(t) = 22 1(t). This 
growth equation has general solution x71(t) = c1e?. The second equation 
then becomes the first order linear equation 


x(t) = 3a4(t) + 3x2(t) 
= 320(t) + 3ce”. 


The integrating factor method applies to find the general solution x2(t) = 
—3c,e7!+ce**. The third and last equation becomes the first order linear 
equation 


x3(t) = 421(t) + 4x9(t) + 423(t) 
= 4x3(t) + 4c,e74 + 4(—3c,e7" + ce**). 


The integrating factor method is repeated to find the general solution 
x3(t) = 4cje”* — 4cge** + cze*. 

In summary, the scalar general solution to the system is given by the 
formulas 


ai(t) = ce”, 

go(t) = —3cye7* + coe, 

x3(t) = 4cje* — 4cge** + cze*. 
Structure of solutions. A system ¥’ = AX for n x n triangular A 
has component solutions x(t), ..., %,(t) given as polynomials times 
exponentials. The exponential factors e@!’, ..., e%"' are expressed in 
terms of the diagonal elements aj1, ..., Gyn of the matrix A. Fewer than 


n distinct exponential factors may appear, due to duplicate diagonal 
elements. These duplications cause the polynomial factors to appear. 
The reader is invited to work out the solution to the system below, 
which has duplicate diagonal entries a1, = ag2 = a33 = 2. 


x1 (t) 2x1(t), 


v(t) = 321(t) + 2ra(t), 
£3(t) 421 (t) + 4x9 (t) + 2x3(t). 
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The solution, given below, has polynomial factors t and t?, appearing 


because of the duplicate diagonal entries 2, 2,2, and only one exponential 


factor e2¢. 


a(t) = ce, 
Bop) (t) => 3c; te” ae ce?" , 
a3(t) = 4cte** + 6cit?e7* + 4cote™ + ce”. 


Conversion to Systems 


Routinely converted to a system of equations of first order are scalar 
second order linear differential equations, systems of scalar second order 
linear differential equations and scalar linear differential equations of 
higher order. 


Scalar second order linear equations. Consider an equation 
au" + bu! + cu = f where a 4 0, b, c, f are allowed to depend on t, 
‘= d/dt. Define the position-velocity substitution 


x(t)=u(t), y(t) =w'(t). 


Then a = uw! = yand y! = wu" = (—bu' —cu+ f)/a = —(b/a)y— (c/a)a + 
f/a. The resulting system is equivalent to the second order equation, in 
the sense that the position-velocity substitution equates solutions of one 
system to the other: 


Te Beh BE 
y@) = ~Don- Oyo +0 


The case of constant coefficients and f a function of t arises often enough 
to isolate the result for further reference. 


Theorem 2 (System Equivalent to Second Order Linear) 
Let a 4 0, b, cbe constants and f(t) continuous. Then aw”’+bu'+cu = f(t) 
is equivalent to the first order system 


aw'(iy= ( Bae Jos ( os ip roe ( res ). 


Converting second order systems to first order systems. A sim- 
ilar position-velocity substitution can be carried out on a system of two 
second order linear differential equations. Assume 


" ! 

ayuy + byuy + cu, = Fis 
" ! 

agus + bou5 +cegu2 = fo. 
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Then the preceding methods for the scalar case give the equivalence 


/ 


ay 0 0 O U1 0 ay 0 0 U1 0 
0 ay 0 O ty _ —Cl —b 0 0 ul ae fi 
0 O ag 0 U2 0 0 0 ag U2 0 
0 0 O ag us 0 0 —C9 —bo us fo 


Coupled spring-mass systems. Springs connecting undamped cou- 
pled masses were considered at the beginning of this chapter, page 754. 
Typical equations are 


myx (t) = —k,x1(t) + ko [x2(t) _ Gilt J 
(1) mgzx9(t) = —ke[wo(t) — vi(t)] + kslag(t) — v2(4)], 
m3x5(t) = —k3[x3(t) cae x2(t)| = kaa3(t). 


The equations can be represented by a second order linear system of 
dimension 3 of the form Mx” = KX, where the position X, the mass 
matrix M and the Hooke’s matrix K are given by the equalities 


LY My 0 0 
x= LQ ; M= 0 meg 0 5 
L3 0 0 M3, 
—(ky + kz) ko 0 
K= ko —(kg + kz) kg 
0 —kg —(k3 + kq) 


Systems of second order linear equations. A second order sys- 
tem MX" = KX + F(t) is called a forced system and F is called the 
external vector force. Such a system can always be converted to a sec- 
ond order system where the mass matrix is the identity, by multiplying 
by M71: 
x" =M'1K¥+M'F(t). 

The benign form ¥” = AX + G(t), where A= M~!K and G = M~!F, 
admits a block matrix conversion into a first order system: 


Ci 4 ae Gee’ x (t) 0 
a (300 ) : (att x) (30 )*( ae): 


Damped second order systems. The addition of a damper to each 
of the masses gives a damped second order system with forcing 


Mx" = BX'+KX+4+F(t). 


In the case of one scalar equation, the matrices M, B, K are constants 
m, —c, —k and the external force is a scalar function f(t), hence the 
system becomes the classical damped spring-mass equation 


ma" + ca’ +kx = f(t). 
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A useful way to write the first order system is to introduce variable 
u = Mx, in order to obtain the two equations 


ii’ = Mx’, i” = Bx’+KxX4+F(t). 


Then a first order system in block matrix form is given by 


Mio \d(xt)\ _ fo|m\(x@ 0 
O|;M J/a\ x(t) )  £#\ K/B z(t) J 7 F(t) } 


The benign form x” = M~! Bx’! + M~!KxX+M~!F (t), obtained by left- 
multiplication by M~', can be similarly written as a first order system 
in block matrix form. 


ay EON Oni 2 X (t) 0 
dt \ x'(t) 7 M"K|M*B x’ (t) ah MF (t) )° 


Higher order linear equations. Every homogeneous nth order 
constant-coefficient linear differential equation 


(n) 


(n—1) 


yo" = Poy +++ + Pn-1Y 


can be converted to a linear homogeneous vector-matrix system 


y 0 1 0 :-- 0 Yy 
F y’ Or Wie al. Agar: y’ 
hs y” a . yf! 
dx ; 
: 0 O O :-- 1 : 
go) Po Pl Po *** Pn-i yr) 


This is a linear system ti’ = Ali where U is the n x 1 column vector 
consisting of y and its successive derivatives, while the n x n matrix A 
is the classical companion matrix of the characteristic polynomial 


2 —1 
r” = po + pir + por* +++++pn—-ir". 


To illustrate, the companion matrix for r+ = a + br + cr? + dr? is 


01 0 0 
00 1 0 
AEN oar) A 
abed 


The preceding companion matrix has the following block matrix form, 
which is representative of all companion matrices. 


_({ 0|/ 1 
A (214-}. 
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Continuous coefficients. It is routinely observed that the methods 
above for conversion to a first order system apply equally as well to 
higher order linear differential equations with continuous coefficients. To 
illustrate, the fourth order linear equation y’” = a(x)y+b(ax)y’ +e(x)y"” + 
d(x)y'” has first order system form ti’ = AU where A is the companion 
matrix for the polynomial r+ = a(x) + b(x)r + c(x)r? + d(x)r3, x held 
fixed. 

Forced higher order linear equations. All that has been said above 
applies equally to a forced linear equation like 


y” = 2y + sin(x)y’ + cos(x)y” + x*y” + f(z). 


It has a conversion to a first order nonhomogeneous linear system 


0 1 O 0 0 y 
email Or 40> ae Oil, <0 42 2 
Gs 0 |? yl! 
2 sing cosx x? f(z) yl" 
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Linear systems. A linear system is a system of differential equa- 
tions of the form 


/ 
Ce SS Oi ve Gintn + fi, 
/ 
Ly = ax ve Gantn + fa, 
(1) : ; ; 
/ 
; a ae vee Qmntn + fm 


where’ = d/dt. Given are the functions a;;(t) and f;(t) on some interval 
a<t<b. The unknowns are the functions x(t), ..., @n(t). 


The system is called homogeneous if all f; = 0, otherwise it is called 
non-homogeneous. 


Matrix Notation for Systems. A non-homogeneous system of 
linear equations (1) is written as the equivalent vector-matrix system 


x’ = A(t)¥ + F(t), 
where 
64 fi Oiy 92> Oa 
x= : . f = : : A= 


In fn Am1l *** Amn 


Existence-uniqueness. The fundamental theorem of Picard and 
Lindeléf applied to the matrix system ¥’ = A(t)X + £(t) says that a 
unique solution x(t) exists for each initial value problem and the solu- 
tion exists on the common interval of continuity of the entries in A(t) 
and f (t). 

Three special results are isolated here, to illustrate how the Picard theory 
is applied to linear systems. 


Theorem 3 (Unique Zero Solution) 
Let A(t) be an m x n matrix with entries continuous on a < t < b. Then 
the initial value problem 


x'= A(t)X, x(0)=0 
has unique solution ¥(t) =0 ona<t <b. 


Theorem 4 (Existence-Uniqueness for Constant Linear Systems) 
Let A(t) = A be an m x n matrix with constant entries and let Xo be any 
m-vector. Then the initial value problem 


/= AX, (0) =X 


has a unique solution X(t) defined for all values of t. 
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Theorem 5 (Uniqueness and Solution Crossings) 

Let A(t) be an m x n matrix with entries continuous on a < t < b and 
assume f (t) is also continuous on a <t <b. If X(t) and ¥(t) are solutions 
of i’ = A(t)ti + F(t) on a < t < b and X(to) = ¥(to) for some fo, 
a <to <b, then X(t) = y(t) fora<t<b. 


Superposition. Linear homogeneous systems have linear structure 
and the solutions to nonhomogeneous systems obey a principle of su- 
perposition. 


Theorem 6 (Linear Structure) 
Let x’ = A(t)X have two solutions X;(t), X2(t). If ky, ko are constants, 


=> 


then X(t) = ky ¥1(t) + ke X(t) is also a solution of x’ = A(t)x. 


The standard basis {w;,}/_,. The Picard-Lindeléf theorem applied 
to initial conditions X(to) = Xo, with ¥o successively set equal to the 
columns of the n x n identity matrix, produces n solutions Ww, ..., 
w,, to the equation ¥’ = A(t)xX, all of which exist on the same interval 
Gt <b: 


The linear structure theorem implies that for any choice of the constants 
C1, +++, Cn, the vector linear combination 


(2) x (t) = cw 1(t) + C2W 2(t) free ft CrW n(t) 


> 


is a solution of x’ = A(t)xX. 
Conversely, if cj, ..., C, are taken to be the components of a given vector 
Xo, then the above linear combination must be the unique solution of 
the initial value problem with xX (to) = Xo. Therefore, all solutions of the 
equation x’ = A(t)X are given by the expression above, where cj, ..., 
Cy are taken to be arbitrary constants. In summary: 


Theorem 7 (Basis) 

The solution set of x’ = A(t)X is an n-dimensional subspace of the vector 
space of all vector-valued functions x(t). Every solution has a unique basis 
expansion (2). 


Theorem 8 (Superposition Principle) 
Let x’ = A(t)¥+£ (t) have a particular solution X(t). If X(t) is any solution 
of X’ = A(t)X¥ +f (t), then X(t) can be decomposed as homogeneous plus 
particular: 

X(t) = Xp(t) + X(t). 
The term X(t) is a certain solution of the homogeneous differential equation 
x’ = A(t)X, which means arbitrary constants cj, c2, ... have been assigned 
certain values. The particular solution X,,(t) can be selected to be free of 
any unresolved or arbitrary constants. 
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Theorem 9 (Difference of Solutions) 
Let x’ = A(t)X +f (t) have two solutions X = u(t) and X = v(t). Define 


= 


y (t) = W(t) — V(t). Then ¥(t) satisfies the homogeneous equation 


General Solution. We explain general solution by example. If a 
formula « = c,e’ + ce" is called a general solution, then it means that 
all possible solutions of the differential equation are expressed by this 
formula. In particular, it means that a given solution can be represented 
by the formula, by specializing values for the constants c,, cg. We expect 
the number of arbitrary constants to be the least possible number. 


The general solution of x’ = A(t) +f (t) is an expression involving arbi- 
trary constants c,, co, ...and certain functions. The expression is often 
given in vector notation, although scalar expressions are commonplace 
and perfectly acceptable. Required is that the expression represents all 
solutions of the differential equation, in the following sense: 


(a) Every assignment of constants produces a solution of 
the differential equation. 


(b) Every possible solution is uniquely obtained from the 
expression by specializing the constants. 


Due to the superposition principle, the constants in the general solution 
are identified as multipliers against solutions of the homogeneous differ- 
ential equation. The general solution has some recognizable structure. 


Theorem 10 (General Solution) 
Let A(t) be n x n and f(t) n x 1, both continuous on an interval a < t < b. 


The linear nonhomogeneous system x’ = A(t)X + f(t) has general solution 
X given by the expression 


¥ = ¥;,(t) + X,(t). 


The term y = Xp,(t) is a general solution of the homogeneous equation 
y’ = A(t)y, in which are to be found n arbitrary constants ci, ..., Cn- 
The term X = X,(t) is a particular solution of x’ = A(t)X +f (t), in which 


there are present no unresolved nor arbitrary constants. 


Recognition of homogeneous solution terms. An expression X 
for the general solution of a nonhomogencous equation X/ = A(t) +£ (t) 
involves arbitrary constants c),..., Cn. It is possible to isolate both terms 
xX, and X, by a simple procedure. 


To find x,, set to zero all arbitrary constants cj, co, ...; the resulting 
expression is free of unresolved and arbitrary constants. 
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To find Xp, we find first the vector solutions y = tiz(t) of y’ = A(t)y, 
which are multiplied by constants cz. Then the general solution x;, of 
the homogeneous equation y’ = A(t)y is given by 


Kye) — cu 1(t) + C2U 2(t) free t Gully t)s 


Use partial derivatives on expression x to find the column vectors 


This technique isolates the vector components of the homogeneous solu- 
tion from any form of the general solution, including scalar formulas for 
the components of x. In any case, the general solution X of the linear 
system ¥/ = A(t)X +f (t) is represented by the expression 


X= cu 1(t) + c2U 2(t) acide CnU p(t) oh X p(t). 


In this expression, each assignment of the constants c1, ..., Cn produces 
a solution of the nonhomogeneous system, and conversely, each possible 
solution of the nonhomogeneous system is obtained by a unique special- 
ization of the constants C1, ..., Cn. 


To illustrate the ideas, consider a 3 x 3 linear system X/ = A(t)X + £ (t) 
with general solution 


given in scalar form by the expressions 


Ci = cet + coe! sibs 
TQ = (cy + ca)et + c3e7", 
g3 = (2co—c)e—* + (4c, — 2c3)e* + 2t. 


To find the vector form of the general solution, we take partial derivatives 


=> 


% OX. 
u;, = —— with respect to the variable names cj, C2, C3: 


Oc 
t -t 
e€ € 0 
uy = ef 5 Uo = et 5 3 = e2t 
Seo det e- Qe2t 
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Finally, 


x= cu 1(t) + c2U 2(t) + c303(t) + Xplt) 


et et 0) t 
= Cj et +c¢9 et + ¢3 e2t + 0 
—e—t + det 9e-* —2e*t Qt 


The expression X = c,U4(t) + cgU2(t) +c3U3(t) +X,(t) satisfies required 
elements (a) and (b) in the definition of general solution. We will develop 
now a way to routinely test the uniqueness requirement in (b). 


Independence. Constants cj, ..., c, in the general solution ¥ = 
Xp, + Xp appear exactly in the expression Xp, which has the form 
Xp = Cy + Cotig +--+ + Cptin. 


A solution ¥ uniquely determines the constants. In particular, the zero 
solution of the homogeneous equation is uniquely represented, which can 
be stated this way: 


cyl, + coGig+---+c,uU,=O0 implies co =c@Q=-:-=c, =0. 


This statement equivalently says that the list of n vector-valued functions 
tii(t),..., t,(t) is linearly independent. 

It is possible to write down a candidate general solution to some 3 x 3 
linear system ¥/ = AX via equations like 


Ly = cet coe! c3e2", 


vq = ce’ + cet 2c3e7", 


x3 = cye’ + ce + 4c3e”. 


This example was constructed to contain a classic mistake, for purposes 
of illustration. 


How can we detect a mistake, given only that this expression is supposed 
to represent the general solution? First of all, we can test that tj, = 
Ox /0c1, U2 = OX /Oc2, U3 = OX /Oc3 are indeed solutions. But to insure 
the unique representation requirement, the vector functions 11, U2, H3 
must be linearly independent. We compute 


et et e2t 
uy = et 5 U2 = et 5 u3 = 967" 
et et Ae?* 


Therefore, WU; = W2, which implies that the functions W1, U2, 3 fail to 
be independent. While is is possible to test. independence by a rudimen- 
tary test based upon the definition, we prefer the following test due to 
Norwegian mathematician N. H. Abel (1802-1829). 
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Theorem 11 (Abel’s Formula and the Wronskian) 

Let X;,(t) = cu1(t) +--+ +cnUpn(t) be a candidate general solution to the 
equation X’ = A(t)xX. In particular, the vector functions Hii (t), ..., H(t) 
are solutions of x’ = A(t)xX. Define the Wronskian by 


w(t) = det((t1(¢)|---|Wn(£))). 
Then Abel’s formula holds: 


w(t) ae, trace(ACs) ids, (45) 8 


In particular, w(t) is either everywhere nonzero or everywhere zero, accord- 
ingly as w(to) £0 or w(to) = 0. 


Theorem 12 (Abel’s Wronskian Test for Independence) 
The vector solutions tij, ..., UW, of X’ = A(t)X are independent if and only 
if the Wronskian w(t) is nonzero for some t = to. 


Clever use of the point to in Abel’s Wronskian test can lead to succinct 
independence tests. For instance, let 


et et e2t 
uy _ ef 5 U2 = ef 5 U3 = Jer 
et et Ae2t 


Then w(t) might appear to be complicated, but w(0) is obviously zero 
because it has two duplicate columns. Therefore, Abel’s Wronskian test 
detects dependence of ti;, 2, tis. 


To illustrate Abel’s Wronskian test when it detects independence, con- 
sider the column vectors 


e e-* 0 
uy = et ; U2 —= et 5 U3 = et 
=p ae Qe" aiDaet 


At t = tp = 0, they become the column vectors 


1 1 0 
gi-{ 1], w=]1], a=] 1 
3 2 —2 


Then w(0) = det((ti1(0)|t2(0)|ti3(0))) = 1 is nonzero, testing indepen- 
dence of Ui, U2, U3. 


>The trace of a square matrix is the sum of its diagonal elements. In literature, 
the formula is called the Abel-Liouville formula. 
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Initial value problems and the rref method. An initial value 
problem is the problem of solving for x, given 
x’ = A(t)¥+£(t), (to) = Xo. 
If a general solution is known, 
X =i (t) +--+ cntin(t) +X,(), 
then the problem of finding X reduces to finding c, ..., Cy, in the relation 
cit; (to) +--+ + Catin(to) + ¥p(to) = Xo. 


This is a matrix equation for the unknown constants c), ..., Cn of the 


= 


form Be = d, where 


B = (tii(to)|---|tin(to)), =] 2 |, d =o — ¥p(to). 
Cn 
The rref-method applies to find ¢. The method is to perform swap, 


combination and multiply operations to C = (Bld) until rref(C) = 
(I|c). 


To illustrate the method, consider the general solution 


Lt = cet + cet +t, 
t= (cy + c2)et + c3e7", 
r3 = (2co—c,)e* + (4c, — 2c3)e7* + 2t. 


We shall solve for ci, €2, cg, given the initial condition 7;(0) = 1, x2(0) 


0, 23(0) = —1. The above relations evaluated at t = 0 give the system 
1 = cye®+coe® + 0, 
= (cy + c2)e° a c3e°, 
—1 = (2c2 —c,)e® + (4c; — 2c3)e° + 0. 


In standard scalar form, this is the 3 x 3 linear system 


Gj Ae . aD => oak 
Ci. os Neos HRS He St), 
3c] i 2c2 an 2c3 = —-l. 


The augmented matrix C’, to be reduced to rref form, is given by 


1 1 
C=|{11 1 0 
3.2 -—2 -I1 
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After the rref process is completed, we obtain 


1 0 0 —-5 
rref(C)=|] 0 10 6 
001 -1 
From this display, we read off the answer cy = —5, co = 6, cz = —1 and 
report the final answer 
xy = —deb+6e* +t, 
gq = eb —e*, 
zz = 17e~t —18e7! + 2t. 


Equilibria. An equilibrium point ¥o of a linear system x’ = A(t)X is 
a constant solution, X(t) = Xo for all t. Mostly, this makes sense when 
A(t) is constant, although the definition applies to continuous systems. 
For a solution ¥ to be constant means X/ = 0, hence all equilibria are 
determined from the equation 


A(t)¥o =O for all ¢. 


This is a homogeneous system of linear algebraic equations to be solved 
for Xq. It is not allowed for the answer Xo to depend on t (if it does, then 
it is not an equilibrium). The theory for a constant matrix A(t) = A 
says that either Xo = 0 is the unique solution or else there are infinitely 
many answers for Xo (the nullity of A is positive). 
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11.4 Matrix Exponential 


The problem 
d ., 


ae 
has a unique solution, according to the Picard-Lindelof theorem. Solve 
the problem n times, when Xg equals a column of the identity matrix, 
and write W1(t), ..., W,(t) for the n solutions so obtained. Define the 
matrix exponential e4! by packaging these n solutions into a matrix: 


(t) = AX(t), (0) = Xo 


et = (Wy(t)|... 


By construction, any possible solution of Lx = AX can be uniquely 


expressed in terms of the matrix exponential eft by the formula 


Matrix Exponential Identities 


Announced here and proved below are various formulas and identities 
for the matrix exponential elt. 


d 
— (eA*) = AcAt Columns satisfy x’ = Ax. 


dt 
° =I Where 0 is the zero matrix. 
Bett — At If AB = BA. 


At Bt _ (A+ B)t 

cAt As _ Alt +s) 
at 

( At)! _ .-At 


eAt = ry (t)Py +--+ 1n(t)Pr 
Ait Aat 
At _ dtp 4 5 (a a 
e€ e€ T= ( eo 
eft — ety + te*(A— yD) 
at a: 
eAt — ett cos ot 1 + © Sle ge al) 
i 3 ante 
n=0 nl 
eAt P-1-Jtp 


lf AB = BA. 


Since At and As commute. 
Equivalently, e4%e~“¢ = J. 


Putzer’s spectral formula — 
see page 784. 


A is 2 x 2, 1 4 2 real. 
A is 2x 2, Ay = Ag real. 


Ais 2x2, Ay=Ao = a+b, 
b> 0. 


Picard series. See page 786. 


Jordan form J = PAP™!. 
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Putzer’s Spectral Formula 


The spectral formula of Putzer applies to a system X¥’ = AX to find its 
general solution. The method uses matrices P,,..., P, constructed from 
A and the eigenvalues \1,...,An of A, matrix multiplication, and the 
solution r(t) of the first order n x n initial value problem 


AD! 0) wos Ws 7 
Le Ae OP ede -'G ; 
BS] OS ee Deiat, ee) |, 
OF 0: yee 


The system is solved by first order scalar methods and back-substitution. 
We will derive the formula separately for the 2 x 2 case (the one used 
most often) and the n x n case. 


Spectral Formula 2 x 2 


The general solution of the 2 x 2 system x’ = Ax is given by the formula 
X(t) = (ri(t)Pi + ra(t) Po) (0), 
where 71, r2, Pi, P. are defined as follows. 
The eigenvalues r = 41, A2 are the two roots of the quadratic equation 
det(A — rI) = 0. 
Define 2 x 2 matrices P,, P2 by the formulas 
Pp=I, Po=A-A\I. 
The functions r1(t), r2(t) are defined by the differential system 


r = M11, PO =, 
r) = oret+ri, r2(0) = 0. 


Proof: The Cayley-Hamilton formula (A — \1/)(A — AzI) = 0 is valid for 
any 2 x 2 matrix A and the two roots r = Aj, A2 of the determinant equality 
det(A—rI) = 0. The Cayley-Hamilton formula is the same as (A—A2)P2 = 0, 
which implies the identity AP: = A2P:. Compute as follows. 


x'(t) = (ri (t)Pi + 73(t)P2) X(0) 
=e 1(t) Py + 11 (t) Po + Agre(t)P2) ¥(0) 
t)A + Agre(t)P2) X (0) 
i(t)A + r2(t)AP2) X (0) 
= An + r2(t)P2) X (0) 
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This proves that X(t) is a solution. Because ®(t) = r1(t)P, + ro(t)P, satisfies 
®(0) = J, then any possible solution of ¥’ = AX can be represented by the 
given formula. The proof is complete. 


Real Distinct Eigenvalues. Suppose A is 2 x 2 having real distinct 
eigenvalues 41, Ag and x(0) is real. Then 


Ait A2T 


Mt are <= E 
At — A2 


and 
Ait _ pAat 


Ay — A2 


X(t) = (ous = (A wt) X (0). 


The matrix exponential formula for real distinct eigenvalues: 


Ait Aat 


—e 
Ay — A2 


At — itp 4 € 


e (A Ail). 


Real Equal Eigenvalues. Suppose A is 2 x 2 having real equal 


eigenvalues Ay = Ay and X(0) is real. Then r; = e*', rg = te!’ and 


rant 


(t) = (eT + te™"(A — ArL)) ¥(0). 
The matrix exponential formula for real equal eigenvalues: 


eAt — Mt 4 te™*(A— AWD). 


Complex Eigenvalues. Suppose A is 2 x 2 having complex eigen- 
values A; = a+ bi with b > 0 and Ap = a — bi. If X(0) is real, then a 
real solution is obtained by taking the real part of the spectral formula. 
This formula is formally identical to the case of real distinct eigenvalues. 
Then 


Ro(X(t)) = (Re(ri(t))E + Re(ra(t)(A — uD) €(0) 
(A =e: ib)1))) £(0) 
(A —al))) 80) 


at Sin bt 


(Re(eloronyr + Re(e 


at Sin bt 


= (ct cos 1 Le 


The matrix exponential formula for complex conjugate eigenvalues: 


in bt 
eft — gat (cos bt I+ _ (A- al))) : 
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How to Remember Putzer’s 2 x 2 Formula. The expressions 


et = 11 (t)I + ro(t)(A— Aid), 


ert Aat 


(1) 


ry(t) =e, ro(t) = 


are enough to generate all three formulas. Fraction rg is the d/dA-Newton 
quotient for r;. It has limit te*!* as \2 > 1, therefore the formula 
includes the case Ay = Az by limiting. If Ay = A2 = a + ib with b > 0, 
then the fraction rz is already real, because it has for z = e*!* and w = d4 


the form 
22. -8ml bt 


a) ee = b 


Taking real parts of expression (1) gives the complex case formula. 


Spectral Formula n x n 
The general solution of ¥’ = AX is given by the formula 


(t) = (ri (t) Pi + ro(t) Po +---+1n(t)Pr) X (0), 


mL 


where rj, 72, .--, Tn, Pi, Po, ..., Py are defined as follows. 


The eigenvalues r = A1,...,An are the roots of the polynomial equation 
det(A — rI) = 0. 
Define n x n matrices P,,..., P, by the formulas 


Pp=I, Py= Pyi(A-Apil) =I (A-AjD), bh =2,...,n. 


The functions r;(t), ..., Tn(t) are defined by the differential system 
ri = AIT 1, T1 (0) = 1, 
TS => Aor + T1; r2(0) = 0, 
fi = Nata eis SO) SO: 


Proof: The Cayley-Hamilton formula (A — A,I)---(A—ApI) = 0 is valid for 
any n x n matrix A and the n roots r = Aj,...,An of the determinant equality 
det(A — rf) = 0. Two facts will be used: (1) The Cayley-Hamilton formula 
implies AP, = A, Pn; (2) The definition of P, implies \;,Py + Pei = AP for 
1<k<n-—1. Compute as follows. 


1] X(t) =(ri()P, +--- +7, (t)P,)X(0) 


2 = ( Akl k (t) Px + S- maf] X(0) 


k=1 k=2 
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n- n-1 
3 = & dere (t)Pe +n (t)AnPn + >> nF] x (0) 
k=1 k=1 
4 = (x rp(t)(AkPe + Presi) +1n (onrs] X (0) 
k=1 
5 = S p(t) AP + rabPa) X (0) 
k=1 
6 =A Sno) X (0) 
k=1 
7 = Ax(t 
Details: |1] Differentiate the formula for X(t). | 2} Use the differential equa- 
tions for r1,...,%. |3 | Split off the last term from the first sum, then re-index 
the last sum. |4| Combine the two sums. |5| Use the recursion for P, and 
the Cayley-Hamilton formula (A — \,J)P, = 0. [6] Factor out A on the left. 


“I 


Apply the definition of X(t). 


This proves that X(t) is a solution. Because ®(t) = )\¢_,rx(t)Py satisfies 
(0) = J, then any possible solution of x’ = AX can be so represented. The 
proof is complete. 


Proofs of Matrix Exponential Properties 


/ 
Verify (e4*) = Ae“*, Let ¥o denote a column of the identity matrix. Define 


X(t) = e4'Xo. Then 


(e4t)'% = X(t) 
= Ax(t) 
= Aetx O- 
Because this identity holds for all columns of the identity matrix, then (e4*)’ and 


Ae“ have identical columns, hence we have proved the identity (Cay = Ae“*, 


Verify AB = BA implies Be“! = eB. Define w(t) = e4*BWo and 
Wo(t) = Bet!Wo. Calculate w{(t) = Awi(t) and w4(t) = BAe4**Wo 
ABe4twW = AWo(t), due to BA = AB. Because W1(0) = W2(0) = Wo, then 
the uniqueness assertion of the Picard-Lindeléf theorem implies that w(t) = 
W(t). Because Wo is any vector, then e4*B = Be4*. The proof is complete. 


Verify eAteBt — elATB)t Tot Xo be a column of the identity matrix. Define 
X(t) = eAteP'Xy and ¥(t) = e(4+)'X9. We must show that X(t) = ¥(t) for 
all t. Define u(t) = e?*X . We will apply the result e4*B = Be4*, valid for 


BA= AB. The details: 
‘(t) 
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We also know that y’(t) = (A+ B)y(t) and since X(0) = ¥ (0) = Xo, then the 
Picard-Lindel6f theorem implies that x(t) = y(t) for all t. This completes the 
proof. 


Verify e4%e4s = eA(t+5) Let ¢ be a variable and consider s fixed. Define 
X(t) = e4%e4°Xq and ¥ (t) = e4¢+)¥q. Then ¥(0) = ¥ (0) and both satisfy the 
differential equation u’(t) = Ati(t). By the uniqueness in the Picard-Lindeléf 
theorem, X(t) = ¥(t), which implies e4%e4* = eA(++s) The proof is complete. 


[o-e) 
: i 
Verify e4! = y Ae The idea of the proof is to apply Picard iteration. 
n! 
n=0 


By definition, the columns of e“? are vector solutions W1(t), ..., Wn(t) whose 
values at t = 0 are the corresponding columns of the n x n identity matrix. 


According to the theory of Picard iterates, a particular iterate is defined by 


t 
Frail) =Fo+ [ AFn(rar 0. 
0 


The vector ¥o equals some column of the identity matrix. The Picard iterates 
can be found explicitly, as follows. 


Vilt) = otf, A¥odr 
= (+ At) ¥o, 
Yo(t) = Yot fs, A¥i(r)dr 


= Yotf, AU + At) Fodr 
I+ At + A?t?/2) ¥o, 


l| 
— 


C4 2 n 1: od 
Yai -= (14 At+ ao +--+ “) Fo. 
The Picard-Lindelof theorem implies that for yg = column k of the identity 


matrix, 


noo 


This being valid for each index k, then the columns of the matrix sum 
N 
a 
m!) 
m=0 
converge as N > co to Wi(t), ..., Wn(t). This implies the matrix identity 
ME he 
— n 
n=0 


The proof is complete. 


Computing e“* 


Theorem 13 (Computing e”’ for J Triangular) 
If J is an upper triangular matrix, then a column u(t) of e7’ can be com- 
puted by solving the system u/(t) = Ju(t), W(0) = V, where V is the 
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corresponding column of the identity matrix. This problem can always be 
solved by first-order scalar methods of growth-decay theory and the inte- 
grating factor method. 


Theorem 14 (Exponential of a Diagonal Matrix) 
For real or complex constants Ay, ..., An, 


ediag(r,....rn)é — diag (Cou Sat et) 


Theorem 15 (Block Diagonal Matrix) 
If A = diag(B,,..., By) and each of By, ..., By is a square matrix, then 


e“! — diag Cay ei Prt) 


Theorem 16 (Complex Exponential) 
Given real a, b, then 


a b : 
—b a an ( cos bt sin 4 
e€ =e . 


—sinbt cos bt 


Exercises 11.4 


Matrix Exponential. #50) = 0 
1 


1. (Picard) Let A be real 2x 2. Write | is triangular so that first-order meth- 
out the two initial value problems | ods can solve the systems. 
which define the columns wj(t), 1 0 
W(t) of e@. 5. A= ( ) 


i: In these exercises A 


2. (Picard) Let A be real 3x3. Write =) 
out the three initial value problems 6. A= ( 0 i: 


0 
0 
which define the columns wj(t), 
W(t), W3(t) of eA. Go nas ( 11 ) 
-A=| 9 9}: 


3. (Definition) Let A be real 2 x 2. 
Show that the solution X(t) =| 8g, 4= ( -iil ). 
e“tiig satisfies ¥/ = AX and 


%(0) = to. Matrix Exponential Identities. 


4. Definition Let A be real n x n. | g, 
Show that the solution x(t) = 
e4txX (0) satisfies X’ = Ax. 10. 

11. 


Matrix Exponential 2 x 2. Find 
et using the formula e4* = (w |W.) | 12- 
and the corresponding systems W4 = | 73. 


Bie th 1 es 3 
AW, Wi(0) = ( 0 iF ws = Awe, | 44. 
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11.5 The Eigenanalysis Method 


The general solution X(t) = e4*X(0) of the linear system 
x’ = AX 


can be obtained entirely by eigenanalysis of the matrix A. A compu- 
tationally useful case is when the n x n matrix A has n independent 
eigenvectors in its list of eigenpairs 


Osta; Domo). seg Aneta) 


It is not required that the eigenvalues 41, ..., An be distinct. The 
eigenvalues can be real or complex. 


The Eigenanalysis Method for a 2 x 2 Matrix 
Suppose that A is 2 x 2 real and has eigenpairs 
(A1,¥1), (2, ¥2), 


with V1, V2 independent. The eigenvalues 41, A2 can be both real. Also, 
they can be a complex conjugate pair Ay = Ag = a+ib with b> 0. 


It will be shown that the general solution of ¥’ = AX can be written as 


X(t) = cee, + C2@*2F 5, 
The details: 
X! = c1(e™")'F 1 + co(€*?")'Vo Differentiate the formula for xX. 
— ce V1 + e902" No¥ 9 
= cet Avy + c9e 2! AV 5 Use AyV 1 = AVq, AoV2q = AVo. 
=A (ae"Wi + ce¥2) Factor A left. 
= Ax Definition of X. 


Let’s rewrite the solution ¥ in the vector-matrix form 


ert Cc 
£(t) = (Wil¥2) ( a ee ( : ). 


Because eigenvectors V1, V2 are assumed independent, then (V1|¥2) is 
invertible and setting t = 0 in the previous display gives 


( = (#1]¥2)1X(0). 


C2 
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Because c1, C2 can be chosen to produce any initial condition x (0), then 
x(t) is the general solution of the system X¥’ = Ax. 


The general solution expressed as ¥(t) = e4*X (0) leads to the exponen- 


tial matrix relation 
Ait 
At 2 tes e 0 yp oye AA 
e = wile ( 0 due ) ta) 


The formula is immediately useful when the eigenpairs are real. 


Complex conjugate eigenvalues. First, eigenpair (Az, V2) is never 
computed or used, because AV, = A,V1 implies AV; = A V1, which 
implies Ap = A; has eigenvector V2 = V1. 


If A is real, then e“* is real, and taking real parts across the formula for 
e4! will give a real formula. Due to the unpleasantness of the complex 
algebra, we will report the answer found, which is real, and then justify 
it with minimal use of complex numbers. 


Define for eigenpair (A1, V1) symbols a, b, P as follows: 
Ay =at+ib, b>0, P= (Re(¥1)|Zm(¥1)). 


Then 
(1) eft — tp ( cos bt sin bt pol 


—sinbt cosbt 


Justification of (1). The formula is established by showing that the matrix 
&(t) on the right satisfies 6(0) = I and ®’ = A®. Then by definition, e4¢ = 
®(t). For exposition, let 


_ pat cosbt sin bt = = 
HU) =e ( —sinbt cos bt ) Pa Oa Cae 


The identity (0) = I verified as follows. 


6(0) = PR(0)P™ 


= Pe (} °) po} 


Write 4; = a+ ib and V1 = Re(V¥i) + iZm(v,). The expansion of eigenpair 
relation AV; = A, ¥V, into real and imaginary parts gives the relation 


A(Re(¥1) + iZm(¥1)) = (a + tb) (Re(¥1) +iZm(¥;)), 


which shows that 
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Then 
@'(t)-'(t) = PR'(t)P7'PR'(t)P 
PR (t)R'(t)P7 


0 b 24 
P at ( | re 


The proof of ®’(t) = A®(t) is complete. 


The formula for e“’ implies that the general solution in this special case 


is 
acy eee 5 5 cos bt sin bt C1 
SO eRe ee) ( —sinbt cos bt ( C2 
The values c;, cg are related to the initial condition x(0) by the matrix 


identity 


C2 


( a = (Re(¥1)|Zm(¥1))71X (0). 


The Eigenanalysis Method for a 3 x 3 Matrix 
Suppose that A is 3 x 3 real and has eigenpairs 
(A1, V1), (A2, V2), (A3, V3); 


with V1, V2, V3 independent. The eigenvalues A;, 2, A3 can be all real. 
Also, there can be one real eigenvalue A3 and a complex conjugate pair 
of eigenvalues Aj = Ag =a+ib with b> 0. 


The general solution of ¥’ = AX can be written as 
xX (t) = ce yy + 02H + C363 ¥ 3. 


The details parallel the 2 x 2 details; they are left as an exercise for the 
reader. 


The solution X¥ is written in vector-matrix form 


emt 9 0 Ci 
X (t) = (¥i|¥2, V3) 0 erat 0 (63) 
0 0 erst C3 


Because the three eigenvectors V1, V2, V3 are assumed independent, 
then (V1|V2|v¥3) is invertible. Setting t = 0 in the previous display gives 


C2 = (¥1,%2,%3) “X(0). 
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Constants c;, C2, c3 can be chosen to produce any initial condition x (0), 
therefore X(t) is the general solution of the 3 x 3 system x’ = AX. There 
is a corresponding exponential matrix relation 


out, DF 0 
et = (Vi|V¥2, V3) 0 er2t 0 (¥1,¥2, 43). 
0 erst 


This formula is normally used when the eigenpairs are real. When there 
is a complex conjugate pair of eigenvalues \1 = Ap = a + ib, b > 0, then 
as was shown in the 2 x 2 case it is possible to extract a real solution 
x from the complex formula and report a real form for the exponential 
matrix: 


e“cosbt e“sinbt 0 


et = P| —e*sinbdt e“cosbt 0 Pol, 
0 0 erst 
P = (Re(¥1)|Zm(V¥ 1), ¥3). 


The Eigenanalysis Method for an n x n Matrix 


The general solution formula and the formula for e4 generalize easily 
from the 2 x 2 and 3 x 3 cases to the general case of an n x n matrix. 


Theorem 17 (The Eigenanalysis Method) 
Let the n x n real matrix A have eigenpairs 


(A1, V1), (A2, V2), sey (Ans Vw) 


with n independent eigenvectors V1, ..., Vn. Then the general solution of 
the linear system x’ = AX is given by 
x)= ev e8" + covge®?? +--+ tonne. 


The vector-matrix form of the general solution is 


Cl 
X(t) = (¥1 ie l¥ n) diag(e™’, scbeeg ent) 
Cn 
This form is real provided all eigenvalues are real. A real form can be 
made from a complex form by following the example of a 3 x 3 matrix 


A. The plan is to list all complex eigenvalues first, in pairs, A, A1, .--, 
Ap; Ap. Then the real eigenvalues r1,...,7¢ are listed, 2p+q =n. Define 


P =(Re(¥1)|Zm(V 1)... | Re(V2p-1)| Zm(V 2p-1)|¥ aptal +++ |¥n); 
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cos bt = sin bt 
—sinbt cos bt 


Ry(t) =e ( 


Then the real vector-matrix form of the general solution is 


). where A+a+ib, b>O0. 


X(t) = Pdiag(Ry,(t),...,Ry,(t),e",...,e) 


and 
A= Pdiag(Ry; O)y.u4 Ri, eee. 


Spectral Theory Methods 


The simplicity of Putzer’s spectral method for computing e“” is appre- 
ciated, but we also recognize that the literature has an algorithm to 
compute e“*, devoid of differential equations, which is of fundamental 
importance in linear algebra. The parallel algorithm computes e4° di- 
rectly from the eigenvalues A; of A and certain products of the nilpotent 
matrices A—.,I. Called spectral formulas, they can be implemented 
in a numerical laboratory or computer algebra system, in order to effi- 
ciently compute e4*, even in the case of multiple eigenvalues. 


Theorem 18 (Computing e“° for Simple Eigenvalues) 


Let the n x n matrix A have n simple eigenvalues Ay, ..., An (possibly 
complex) and define constant matrices Qi, ee Qn, by the formulas 
= A-xXI 


Then 
eft = evga, Aeieany the ertQ,,. 


Theorem 19 (Computing e“’ for Multiple Eigenvalues) 

Let the n x n matrix A have k distinct eigenvalues A;, ..., Ax of algebraic 
multiplicities m1, ..., Mz. Let p(A) = det(A — AJ) and define polynomials 
ay(A), ..-, @e(A) by the partial fraction identity 


1 aw(\) eer ay, (A) 
pA) (A= A1)™ (A= Ag) 


Define constant matrices Qi, ahs QO: by the formulas 


Q;=a,;(AIiz,(A-Ad)™, g=1)...,k: 
Then 
p 


(2) a se aS A-XI)I 5. 
i=1 j=0 
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Proof: Let N; = Q;(A —XI),1<i<k. We first prove 


Lemma 1 (Properties) 


~-Qit---+Qr= 

2. Q:Q;i =Qi, 

3. Q;Q,; =0 fori Fj, 
4. NiN; =0 fori 4 j, 
5. N™ —0, 

6. A= (AQ; + Ni). 


The proof of 1 follows from clearing fractions in the partial fraction expansion 
of 1/p(A): 


The projection property 2 follows by multiplication of identity 1 by Q; and 
then using 2. 

The proof of 3 starts by observing that Q; and Q ; together contain all the fac- 
tors of p(A), therefore Q:Q; = q(A)p(A) for some polynomial g. The Cayley- 
Hamilton theorem p(A) = 0 finishes the proof. 

To prove 4, write NiN; =(A-A,I)(A- \,1)Q:Q,; and apply 3. 

m: — Q, (from 2) to write N™ = (A—A,I)™Q, = p(A) = 9. 
To prove 6, multiply 1 by A and rearrange as follows: 


A ei AQ: , 
5 NiQi TT (A 7 AL)Q;: 


= Yer rGs+N; 


To prove 5, use Q 


To prove (2), multiply 1 by e4* and compute as follows: 
$50. 
— mee Q je" 
_ yk Q eriltt+(A-Ail)t 
au 
= ae Q ertelA-ADt 
_ ee Qjcrite@i(A-At 
a sk QjecriteNit 
= Th Qie™* NI (A= AD) 


j! 


Solving Planar Systems x’(t) = Ax (t) 


A 2x 2 real system x/(t) = AX(t) can be solved in terms of the roots of 
the characteristic equation det(A — AJ) = 0 and the real matrix A. 


Theorem 20 (Planar System, Putzer’s Spectral Formula) 
Consider the real planar system x/(t) = AxX(t). Let Ay, A2 be the roots of 
the characteristic equation det(A— AI) = 0. The real general solution x (t) 


is given by the formula 
X(t) = e4'X (0) 
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where the 2 x 2 exponential matrix e4* is given as follows. 
Agt _ Ait 
Real A; # Ao gibt eT) 
A2 - AL 
Real Ay = Az eft — et'T + te**(A — ALD). 
_ at a: bt 
Complex A; = Ao, eft — e% cos bt I + eT —al). 


Ay =at+bi,b>0 


Proof: The formulas are from Putzer’s algorithm, or equivalently, from the 
spectral formulas, with rearranged terms. The complex case is formally the 
real part of the distinct root case when \2 = A,. The spectral formula is the 
analog of the second order equation formulas, Theorem 45 in Chapter 5. 


Illustrations. Typical cases are represented by the following 2 x 2 
matrices A, which correspond to roots 1, A2 of the characteristic equa- 
tion det(A — AI) = 0 which are real distinct, real double or complex 


conjugate. The solution ¥(t) = e4*X(0) is given here in two forms, by 


writing e“* using |1] a spectral formula and |2) Putzer’s spectral 
formula. 
Ay = 5, Ag = 2 Real distinct roots. 


56 or oe 
A= =P 1) a2 2 iis as 
-—6 8 3 \-6 6 -3\-6 3 


Qt 5t 
At bt e~" —e€ —6 3 
2) e" =e ar re & ? 


Ay = A2 =3 Real double root. 


2 1 At _ ,3t 1, AE 
a= (4 :) 1|e*“ =e ae ')) 


1 1 
Ay = Ag = 243i Complex conjugate roots. 
a=(3 3] He =ore (> (5 a] 
2) e4! = ec cos 3tI + ea & 


The complex example is typical for real n x n matrices A with a complex 
conjugate pair of eigenvalues \y = A». Then Q> = Q,. The result is 
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that A2 is not used and we write instead a simpler expression using the 
college algebra equality z+ 7 = 2 Re(z): 
evlg, +e'Qs =2Re (e"Q 1) . 


This observation explains why e“¢ is real when A is real, by pairing 
complex conjugate eigenvalues in the spectral formula. 
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11.6 Jordan Form and Ejigenanalysis 


Generalized Eigenanalysis 


The main result is Jordan’s decomposition 
A=PJP™, 


valid for any real or complex square matrix A. We describe here how 
to compute the invertible matrix P of generalized eigenvectors and the 
upper triangular matrix J, called a Jordan form of A. 


Jordan block. Anmxm upper triangular matrix B(\, m) is called a 
Jordan block provided all m diagonal elements are the same eigenvalue 
X and all super-diagonal elements are one: 


A 1 0 -:- 0 0 
OA 1: 0 0 
BOS) = a a (m x m matrix) 
0 0 0 rA 1 
0 0 0 0 A 


Jordan form. Given an n X n matrix A, a Jordan form J for A is 
a block diagonal matrix 


J= diag(B(A1, m1), BU, mz), Pa » BAR, Me)), 


where Ay, ..., Ax are eigenvalues of A (duplicates possible) and m1 + 
--> +m, = n. Because the eigenvalues of A are on the diagonal of J, 
then A has exactly k eigenpairs. If k <n, then A is non-diagonalizable. 
The relation A = PJP! is called a Jordan decomposition of A. 
Invertible matrix P is called the matrix of generalized eigenvectors 
of A. It defines a coordinate system X = Py in which the vector function 
xX — AX is transformed to the simpler vector function y > Jy. 


If equal eigenvalues are adjacent in J, then Jordan blocks with equal 
diagonal entries will be adjacent. Zeros can appear on the super-diagonal 
of J, because adjacent Jordan blocks join on the super-diagonal with a 
zero. A complete specification of how to build J from A appears below. 


Decoding a Jordan Decomposition A = PJP™'. If Jisasingle 
Jordan block, J = B(A,m), then P = (¥1|...|¥m) and AP = PJ means 


Ay. 3S OH 
AV2 = AV24+V1, 


AVm = AVm+Vm-1- 
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The exploded view of the relation AP = PB(A,m) is called a Jordan 
chain. The formulas can be compacted via matrix N = A— AI into the 
recursion 

N¥,=0, N¥2.=¥1,...,N¥m=V¥m-1- 


The first vector V1 is an eigenvector. The remaining vectors V2, ..., Vim 
are not eigenvectors, they are called generalized eigenvectors. A 
similar formula can be written for each distinct eigenvalue of a matrix A. 
The collection of formulas are called Jordan chain relations. A given 
eigenvalue may appear multiple times in the chain relations, due to the 
appearance of two or more Jordan blocks with the same eigenvalue. 


Theorem 21 (Jordan Decomposition) 
Every n x n matrix A has a Jordan decomposition A = PJP-!. 


Proof: The result holds by default for 1 x 1 matrices. Assume the result holds 
for all k x k matrices, k <n. The proof proceeds by induction on n. 


The induction assumes that for any k x k matrix A, there is a Jordan decom- 
position A = P.JP~'. Then the columns of P satisfy Jordan chain relations 


oj oj ae ; ot =I 
AX} =AXP+ KI, Gg >1, AR; = riX;- 


Conversely, if the Jordan chain relations are satisfied for k independent vectors 
{x/}, then the vectors form the columns of an invertible matrix P such that 
A= PJP with J in Jordan form. The induction step centers upon producing 
the chain relations and proving that the n vectors are independent. 


Let B be n x n and Ao an eigenvalue of B. The Jordan chain relations hold for 
A= B if and only if they hold for A = B— Aol. Without loss of generality, we 
can assume 0 is an eigenvalue of B. 


Because B has 0 as an eigenvalue, then p = dim(kernel(B)) > 0 and k = 
dim(Image(B)) < n, with p+k=n. Ifk =0, then B = 0, which is a Jordan 
form, and there is nothing to prove. Assume henceforth p and k positive. Let 
S = (col(B,i1)|...|col(B,ix)) denote the matrix of pivot columns 1j,... ,7% 
of B. The pivot columns are known to span Image(B). Let A be the k x k 
basis representation matrix defined by the equation BS = S'A, or equivalently, 
Bcol(S,j) = uae a;; col(S,i). The induction hypothesis applied to A implies 
there is a basis of k-vectors satisfying Jordan chain relations 


oj Sod Ge : = = 
AX? =A;X7+xX77, j>1, Axt = A;x). 


The values \;, i = 1,...,p, are the distinct eigenvalues of A. Apply S to these 

equations to obtain for the n-vectors y} = Sx? the Jordan chain relations 
BY} =M¥{+ HP", G>1, BY; =X. 

Because S has independent columns and the k-vectors x! are independent, then 

the n-vectors y7 are independent. 

The plan is to isolate the chains for eigenvalue zero, then extend these chains 


by one vector. Then 1-chains will be constructed from eigenpairs for eigenvalue 
zero to make n generalized eigenvectors. 
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Suppose gq values of 7 satisfy 4; = 0. We allow q = 0. For simplicity, assume such 
values i are i = 1,...,q. The key formula y? = SxX/ implies y? is in Image(B), 
while BY; = Ax¥; implies yj,...,¥{ are in kernel(B). Each eigenvector y} 
starts a Jordan chain ending in yn), Then® the equation BU = yo has 
an n-vector solution Gd. We label i = york, Because A; = 0, then BU = 


Aya + ym results in an extended Jordan chain 


BY; = MY; 

BY? = MY} = ve 
By i yym® 4 gm@- 
By = hy ge ee 


Let’s extend the independent set {y}}7_, to a basis of kernel(B) by adding 
s = n—k—q additional independent vectors V1,...,V5. This basis consists 
of eigenvectors of B for eigenvalue 0. Then the set of n vectors V;, ¥? for 
l<r<s,l<i<p,1<j < m(t)+1 consists of eigenvectors of B and vectors 
that satisfy Jordan chain relations. These vectors are columns of a matrix P 
that satisfies BP = PI where 7 is a Jordan form. 


To prove P invertible, assume a linear combination of the columns of P is zero: 


p mi) sg m(i41 ig 
dD ut DO ay + avi =0. 
i=q41 j=l i=1 j=l i=l 


Apply B to this equation. Because BW =0 for any W in kernel(B), then 


po mi) gmt 
So So uBy +S" SO By] =06. 
i=q+1 j=1 i=1 j=2 


The Jordan chain relations imply that the k vectors By! in the linear com- 
bination consist of ay! + ys MY}, 4=@4+1,...,p, 7 = 2,...,m(4), plus 
the vectors ¥!, 1<i<q,1< j< m(t). Independence of the original k vec- 
tors {yi} plus A; 4 0 for i > q implies this new set is independent. Then all 
coefficients in the linear combination are zero. 

The first linear combination then reduces to 74_, bl¥!} + 338_, av; = 0. In- 
dependence of the constructed basis for kernel(B) implies b} = 0 for 1 < i < q 
and c; = 0 for 1 <i<_s. Therefore, the columns of P are independent. The 
induction is complete. 


Geometric and algebraic multiplicity. The geometric multi- 
plicity is defined by GeoMult(A) = dim(kernel(A — AJ)), which is the 
number of basis vectors in a solution to (A—AJ)X = 0, or, equivalently, 
the number of free variables. The algebraic multiplicity is the integer 
k = AlgMult()) such that (r — A)* divides the characteristic polynomial 
det(A — AI), but larger powers do not. 


The n-vector UH is constructed by setting G@ = 6, then copy components of k-vector 
>m(i) om(i) - 


X;" into pivot locations: row(U,i;) = row(X;"°",j), 9 =1,...,k. 
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Theorem 22 (Algebraic and Geometric Multiplicity) 
Let A be a square real or complex matrix. Then 


(1) 1 < GeoMult(\) < AlgMult()). 


In addition, there are the following relationships between the Jordan form J 
and algebraic and geometric multiplicities. 


GeoMult()) Equals the number of Jordan blocks in J with eigen- 
value X, 

AlgMult(\) Equals the number of times A is repeated along the 
diagonal of J. 


Proof: Let d = GeoMult(\,). Construct a basis v1,...,U, of R” such that 
V1,...,Ua is a basis for kernel(A — Aol). Define S = (vi|...|vn) and B = 


S-!AS. The first d columns of AS are \ov1, ..., Apvg. Then B = ( Aot “ ) 


for some matrices C and D. Cofactor expansion implies some polynomial g 
satisfies 
det(A — AI) = det($(B — XI)S~1) = det(B — XT) = (A— Xo) %g(A) 


and therefore d < AlgMult(Ao). Other details of proof are left to the reader. 


Chains of generalized eigenvectors. Given an eigenvalue \ of 
the matrix A, the topic of generalized eigenanalysis determines a Jordan 
block B(A,m) in J by finding an m-chain of generalized eigenvectors 


V1, ...; Vm, which appear as columns of P in the relation A = PJP™!. 

The very first vector V, of the chain is an eigenvector, (A — AJ)v; = 0. 

The others Vo, ..., Vz are not eigenvectors but satisfy 
(A-AlDVo=Vi1, .-. » (A-ADVm=Vm-1.- 


Implied by the term m-chain is insolvability of (A — AI)X = Vm. The 
chain size m is subject to the inequality 1 << m < AlgMult()). 

The Jordan form J may contain several Jordan blocks for one eigenvalue 
d. To illustrate, if J has only one eigenvalue \ and AlgMult(\) = 3, 
then J might be constructed as follows: 


A 0 0 

J =diag(B(A,1),B(,1),BX,1)) = | 0 A 0 |, 
0 0A 
A 0 0 

J = diag(B(A, 1), B(A,2)) S Wore at 3 
0 0A 
A 1 0 

J = B(X,3) =!10)A1 
0 0A 


The three generalized eigenvectors for this example correspond to 
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+ Three 1-chains, 


+ One 1-chain and one 2-chain, 


4+ One 3-chain. 


= 
I 
coOoYoOoOY-oo»y 
OyvyrFOYoTOYO 
~rRPOoOY,rFOYoO 


Computing m-chains. Let us fix the discussion to an eigenvalue 
of A. Define N = A— XI and p = AlgMult()). 


To compute an m-chain, start with an eigenvector V; and solve recur- 
sively by rref methods Nvj+41 = Vv; until there fails to be a solution. 
This must seemingly be done for all possible choices of ¥;! The search for 
m-chains terminates when p independent generalized eigenvectors have 
been calculated. 


If A has an essentially unique eigenpair (A, ¥1), then this process termi- 
nates immediately with an m-chain where m = p. The chain produces 
one Jordan block B(A,m) and the generalized eigenvectors V1, ..., Vim 
are recorded into the matrix P. 


If Uy , U2 form a basis for the eigenvectors of A corresponding to A, then 
the problem NX = 0 has 2 free variables. Therefore, we seek to find an 
my,-chain and an m2-chain such that m +m, = p, corresponding to two 
Jordan blocks B(A,m 1) and B(A, m2). 


To understand the logic applied here, the reader should verify that for 
N = diag(B(0,m 1), B(0,mz2),..., B(0,mp)) the problem NX = 0 has 
k free variables, because NV is already in rref form. These remarks 
imply that a k-dimensional basis of eigenvectors of A for eigenvalue 
causes a search for m;-chains, 1 <i < k, such that m; +--- +m, = p, 
corresponding to k Jordan blocks B(A,m 4), ..., B(A, mx). 


A common naive approach for computing generalized eigenvectors can 
be illustrated by letting 


| te i 
A={010]|, w=] -1 |], w=| 1 
001 1 =i 


Matrix A has one eigenvalue \ = 1 and two eigenpairs (1, t;), (1, ti2). 
Starting a chain calculation with V1 equal to either GW, or U2 gives a 
1-chain. This naive approach leads to only two independent generalized 
eigenvectors. However, the calculation must proceed until three inde- 
pendent generalized eigenvectors have been computed. To resolve the 
trouble, keep a 1-chain, say the one generated by 1, and start a new 
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chain calculation using V; = ajU 1+ adgti2. Adjust the values of a4, ag 
until a 2-chain has been computed: 


1 ay 0 
—a, + ao ~ 0 
a, — ag 0 


0 1 ay 
(A-Al|¥i)=] 0 0 
0 0 


SO: 


1 
0 0 0 |, 
0 0 0 
provided a, — ag = 0. Choose a; = ag = 1 to make V; = U, + tin # 0 
and solve for V2 = (0,1,0 ): Then Uy, is a l-chain and Vj, Vo isa 
2-chain. The generalized eigenvectors U1, V1, V2 are independent and 
form the columns of P while J = diag(B(A, 1), B(A,2)) (recall A = 1). 
We justify A = PJP! by testing AP = PJ, using the formulas 


AO 30 110 
J=}0.2A1]|, P=[{-101 
O20 100 


Jordan Decomposition using maple 


Displayed here is maple code which applied to the matrix 


4 -2 5 
A=]|-2 4 -8 
0 O 2 


produces the Jordan decomposition 


; je ee 6.00) af Bis e 16 
ASTI = es | @. 2571 7 a <>. os 
00 4 002 0 Oo 4 


A d= Matrix¢[[4, <2; 5]. [-2) 4, -31,. 10; 0, 213 

factor (LinearAlgebra[CharacteristicPolynomial] (A, lambda) ) ; 
# Answer == (lambda-6)*(lambda-2) ~2 
J,P:=LinearAlgebra[JordanForm] (A,output=[’J’,’Q’]); 
zero:=A.P-P.J; # zero matrix expected 


Number of Jordan Blocks 


In calculating generalized eigenvectors of A for eigenvalue A, it is pos- 
sible to decide in advance how many Jordan chains of size k should be 
computed. A practical consequence is to organize the computation for 
certain chain sizes. 


Theorem 23 (Number of Jordan Blocks) 

Given eigenvalue of A, define N = A — XI, k(j) = dim(kernel(N’)). 
Let p be the least integer such that N? = N?+!. Then the Jordan form of 
A has 2k(j — 1) —k(j — 2) —k(j) Jordan blocks B(A, 7 — 1), 7 =3,...,p. 
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The proof of the theorem is in the exercises, where more detail appears 
for p = 1 and p= 2. Complete results are in the maple code below. 


An Illustration. This example is a 5x5 matrix A with one eigenvalue 
A = 2 of multiplicity 5. Let s(j7) = number of j x 7 Jordan blocks. 


3-1 1 0 0 A el oil 0 0 

2 O01 1 0 De Den li 1 0 

A= 1 -1 2 10], N=A-2l= 1 -l1 0 1 0 
—1 10 2 1 —1 10 O01 

-3 3 0 -2 3 -3 3 0 -2 1 


Then N? = N+ = N° = 0 implies k(3) = k(4) = k(5) = 5. Further, 
k(2) = 4, k(1) = 2. Then s(5) = s(4) = s(2).=1, 60) =0, 
which implies one block of each size 2 and 3. 


= 
wW 
— 
w 
wn 
I 


Some maple code automates the investigation: 


with(LinearAlgebra) : 

A := Matrix([ 

[35 Rag he 05 Ole Bi Os Des 1g Os 
[ty Sd Bel Oly ede Oy. Be ST 


[HS 85 Oy HD 3): 8 

lambda:=2; 
n:=RowDimension(A) ;N:=A-lLambda*IdentityMatrix(n) ; 
for j from 1 to n do 

k[j] :=n-Rank(N*j); od: 

for p from n to 2 by -1 do 

if (k[p]<>k[p-1])then break; fi: od; 
txt:=(j,x)->printf (‘if ‘ (x=1, 
cat("B(lambda,",j,") occurs 1 time\n"), 
cat("B(lambda,",j,") occurs ",x," times\n"))): 
printf ("lambda=%d, nilpotency=%d\n",, lambda, p) ; 
if(p=1) then txt(1,k[1]); else 

txt (p,k[p]-k[p-1]); 

for j from p to 3 by -1 do 

txt (j-1,2*k[j-1]-k[j-2]-k[j]): od: 

txt (1,2*k[1]-k[2]); 
fi: 
#lambda=2, nilpotency=3 
#B(lambda,3) occurs 1 time 
#B(lambda,2) occurs 1 time 
#B(lambda,1) occurs 0 times 
J,P:=JordanForm(A,output=[’J’,’Q’])}: 
# Answer check for the maple code 
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210 00 G i “Ook wo 
i a ee ,|-4 2 2-2 2 
P5068. Or We PSS ee a et 
00021 SN areas A 245. 
00002 ae ae ae oe 


Numerical Instability 


1 1 


The matrix A = ( a has two possible Jordan forms 


1 1 
tat e=0, 


fee 4 e>0. 


J(e) = 


Os ewe 


When ¢ & 0, then numerical algorithms become unstable, unable to lock 
onto the correct Jordan form. Briefly, limz_.9 J(e) € J(0). 


The Real Jordan Form of A 


Given a real matrix A, generalized eigenanalysis seeks to find a real 
invertible matrix P and a real upper triangular block matrix R such 
that A= PRP7!. 


If \ is a real eigenvalue of A, then a real Jordan block is a matrix 


0410 -:- O 
B=diag(\,...,A)+N, N= Be. Sele oct 
000 .--- J 
0 0 0 0 


If \=a+ib is a complex eigenvalue of A, then symbols A, 1 and 0 are 
a b ‘ 

Sie , Z = diag(1, 1) 

and O = diag(0,0). The corresponding 2m x2m real Jordan block matrix 

is given by the formula 


replaced respectively by 2 x 2 real matrices A = 


B=diag(A,...,A)+N, N= 


SS 
CG <2 
OQ: 
SS 
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Direct Sum Decomposition 


The generalized eigenspace of eigenvalue \ of an n xn matrix A is the 
subspace kernel((A — AJ)?) where p = AlgMult(\). We state without 
proof the main result for generalized eigenspace bases, because details 
appear in the exercises. A proof is included for the direct sum decompo- 
sition, even though Putzer’s spectral theory independently produces the 
same decomposition. 


Theorem 24 (Generalized Eigenspace Basis) 

The subspace kernel((A — AZ)*), k = AlgMult(\) has a k-dimensional 
basis whose vectors are the columns of P corresponding to blocks B(A, 7) 
of J, in Jordan decomposition A = PJP7!. 


Theorem 25 (Direct Sum Decomposition) 
Given nxn matrix A and distinct eigenvalues Ay, ..., Ax, m1 = AlgMult(,;), 
.., Ng = AlgMult(\;), then A induces a direct sum decomposition 


C” = kernel((A — 1J)™ @---@ kernel((A — Ax J)”. 


This equation means that each complex vector X in C” can be uniquely 
written as 


X=X,+--- +X, 


where each x; belongs to kernel ((A — \;)™), i =1,...,k. 


Proof: The previous theorem implies there is a basis of dimension n; for E; = 
kernel((A — A;J)"'), i=1,...,k. Because ny +---+n, =n, then there are n 
vectors in the union of these bases. The independence test for these n vectors 
amounts to showing that ¥; +---+X, =0 with X; in Ej, i =1,...,k, implies 
all X; =0. This will be true provided EB; Ej = {0} for i Aj. 

Let’s assume a Jordan decomposition A = PJ P~'. If ¥ is common to both F; 
and £;, then basis expansion of X in both subspaces implies a linear combina- 
tion of the columns of P is zero, which by independence of the columns of P 
implies ¥ = 0. 


The proof is complete. 


Computing Exponential Matrices 

Discussed here are methods for finding a real exponential matrix e4¢ 
when A is real. Two formulas are given, one for a real eigenvalue and 
one for a complex eigenvalue. These formulas supplement the spectral 
formulas given earlier in the text. 
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Nilpotent matrices. A matrix N which satisfies N? = 0 for some 
integer p is called nilpotent. The least integer p for which N? = 0 is 
called the nilpotency of N. A nilpotent matrix N has a finite exponen- 
tial series: 


= 
(p— 1)! 


t2 
eNi_T4Nt4 We tpt de et 


If N = B(A, p) — AT, then the finite sum has a splendidly simple expres- 
sion. Due to e+N* = eeN* this proves the following result. 
Theorem 26 (Exponential of a Jordan Block Matrix) 

If A is real and 


A 1 0 0 0 
B= obo (p X p matrix) 
0 0 0 Xl 
O00 0 A 
then 
42 ¢Pp—2 tp-1 
healt Sy (2! Hw! 
is end ee a 
OD} wOy ares 1 t 
0 0 0 0 1 


The equality also holds if X is a complex number, in which case both sides 
of the equation are complex. 


Real Exponentials for Complex A. A Jordan decomposition 
A = PJP~—!, in which A has only real eigenvalues, has real general- 
ized eigenvectors appearing as columns in the matrix P, in the natural 
order given in J. When A = a+ ib is complex, b > 0, then the real 
and imaginary parts of each generalized eigenvector are entered pairwise 
into P; the conjugate eigenvalue \ = a — ib is skipped. The complex 
entry along the diagonal of J is changed into a 2 x 2 matrix under the 
correspondence 


The result is a real matrix P and a real block upper triangular matrix J 
which satisfy A = PJP7!. 


Theorem 27 (Real Block Diagonal Matrix, Eigenvalue a + ib) 
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a b 


Let A= he wel T = diag(1,1) and O = diag(0,0). Consider a real 
Jordan block matrix of dimension 2m x 2m given by the formula 
A ZO .::-. O O 
a ea ae ae 
O O O A fF 
O O O O A 


cos bt sin bt 
eS ( —sinbt cos bt ) men 


2 m—2 m—-1 
RR GR Gay may 
Prego a ee : 
OOO Oo. Se Kk tR 
OP OP. SO) ate O R 


Solving x’ = AX. The solution X(t) = e“*X(0) must be real if A is 


real. The real solution can be expressed as X(t) = Py (t) where y’(t) = 
Ry (t) and R is a real Jordan form of A, containing real Jordan blocks 
Bi, ..., By, down its diagonal. Theorems above provide explicit formulas 


for the block matrices e?* in the relation 


e® — diag ee ae est) 
The resulting formula 
X(t) = Pe*P-1¥ (0) 


contains only real numbers, real exponentials, plus sine and cosine terms, 
which are possibly multiplied by polynomials in t. 
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Exercises 11.6 


Jordan block. Write out explicitly. 
1. 
2. 
3. 


A. 


Jordan form. Which are Jordan 
forms and which are not? Explain. 


5. 
6. 
7. 


8. 


Decoding A = PJP™'. Decode 
A = PJP in each case, displaying 
explicitly the Jordan chain relations. 


9. 
10. 
11. 


12. 


Deter- 
multiplicity 


Geometric multiplicity. 
mine the geometric 


GeoMult()). 
13. 
14. 
15. 


16. 


Algebraic multiplicity. Determine 
the algebraic multiplicity AlgMult()). 


17. 
18. 


19. 
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20. 


Generalized eigenvectors. Find all 
generalized eigenvectors and represent 
Aa PPS 


21. 
22. 


32. 


Computing m-chains. Find the Jor- 
dan chains for the given eigenvalue. 


33. 


40. 


Jordan Decomposition. Use maple 
to find the Jordan decomposition. 


41. 
42. 
43. 


11.6 Jordan Form and Ejigenanalysis 


44. 
45. 
46. 
AZ. 


48. 


Number of Jordan Blocks. Outlined 
here is the derivation of 


Definitions: 


e s(j)= number of blocks B(A, 7) 
e N=A-AXI 
e k(j) = dim(kernel(V’)) 


e L; =kernel(N/~')+ relative to 
kernel(N/) 


© (j) = dim(L;) 
e p minimizes 


kernel(N?) = kernel(N?*") 


A9. Verify k(j) < k(j +1) from 


kernel(N’) Cc kernel(N’*"). 


50. Verify the direct sum formula 


kernel(N’) = kernel(N!~') @ L;. 
Then k(j) = k(j — 1) + €(). 


Given Ni¥ = 0, NI-!¥ 4 O, 

define V; = NJ-'¥, i = 1,...,). 
Show that these are independent 
vectors satisfying Jordan chain re- 
lations NV = 0, NV isi = Vi- 


dl. 


52. A block B(A,p) corresponds to 
a Jordan chain Vj, ..., Wp con- 
structed from the Jordan decom- 
position. Use NJ-!¥; = Vi 
and kernel(N?) = kernel(N?*!) 
to show that the number of such 
blocks B(A,p) is &(p). Then for 


p> 1, s(p) = k(p) — k(p— 1). 
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53. Show that ¢(j — 1) — &(j) is the 
number of blocks B(A,j) for 2 < 
j <p. Then 


s(j) = 2k — 1) 


k(j) — kG — 2). 


54. Test the formulas above on the 
special matrices 
A = diag(B(A,1), BOA, 1), BO, 1)), 
A = diag(B(A, 1), BOA, 2), BO, 3)), 
A = diag(B(A, 1), BOA, 3), BO, 3)), 
A = diag(B(A, 1), BOA, 1), BO, 3)), 


Generalized Eigenspace Basis. 

Let A be n x n with distinct eigenval- 
ues Nis y= AlgMult();) and Ei = 
kernel((A — \;J)"’), 1 =1,...,k. As- 
sume a Jordan decomposition A = 
PJP-. 


55. Let Jordan block B(A, j) appear in 
J. Prove that a Jordan chain cor- 
responding to this block is a set of 
j independent columns of P. 


. Let By be the union of all 
columns of P originating from Jor- 
dan chains associated with Jordan 
blocks B(A, 7). Prove that B) is an 
independent set. 


57. Verify that By has AlgMult(\) ba- 


sis elements. 


58. Prove that E; = span(B,,) and 


Numerical Instability. Show directly 
that lim.4o J(e) # J(0). 


59. 


62. 


Direct Sum Decomposition. Dis- 
play the direct sum decomposition. 


63. 
64. 
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65. 
66. 
67. 
68. 
69. 


70. 


Exponential Matrices. Compute the 
exponential matrix on paper and then 
check the answer using maple. 


71. 
72. 
73. 
74. 
75. 
76. 
77. 


78. 


Nilpotent matrices. Find the nilpo- 
tency of N. 


79. 


80. 
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81. 


82. 


Real Exponentials. Compute the 
real exponential e4* on paper. Check 
the answer in maple. 


83. 
84. 
85. 


86. 


Real Jordan Form. Find the real Jor- 
dan form. 


87. 
88. 
89. 


90. 


Solving x’ = AX. Solve the differ- 
ential equation. 


91. 
92. 
93. 


94. 
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11.7 Nonhomogeneous Linear Systems 


Variation of Parameters 


The method of variation of parameters is a general method for 
solving a linear nonhomogeneous system 


x’ = AX + F(t). 


Historically, it was a trial solution method, whereby the nonhomogeneous 
system is solved using a trial solution of the form 


¥ (t) = e4* ¥o(t). 


In this formula, ¥o(t) is a vector function to be determined. The method 
is imagined to originate by varying Xo in the general solution x(t) = 
e4' ¥9 of the linear homogenous system ¥! = AX. Hence was coined the 
names variation of parameters and variation of constants. 


Modern use of variation of parameters is through a formula, memorized 
for routine use. 


Theorem 28 (Variation of Parameters for Systems) 
Let A be a constant n x nm matrix and F(t) a continuous function near 
t = to. The unique solution x(t) of the matrix initial value problem 


X(t) = AX(t) + F(t), (to) = Xo, 


is given by the variation of parameters formula 


t oa 
(1) X(t) = e4MtXy + et | TAB (r)dr. 
to 


Proof of (1). Define 
t > 
u(t) — Xo +/ e TAR (r)dr. 
to 


To show (1) holds, we must verify X(t) = e4‘W(t). First, the function W(t) is 
differentiable with continuous derivative e~‘4F (t), by the fundamental theorem 
of calculus applied to each of its components. The product rule of calculus 
applies to give 
X(t) = (et) H(t) + e4*a(t) 
= Ae*tii(t) + e4te-4*F (t) 
) 


Therefore, X(t) satisfies the differential equation X/ = AX + F(t). Because 
U (to) = Xo, then X(to) = Xo, which shows the initial condition is also satisfied. 
The proof is complete. 
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Undetermined Coefficients 


The trial solution method known as the method of undetermined coef- 
ficients can be applied to vector-matrix systems X’ = AX + F(t) when 
the components of F are sums of terms of the form 


(polynomial in t)e“‘(cos(bt) or sin(bt)). 


Such terms are known as Euler solution atoms. It is usually efficient 
to write F in terms of the columns 6), ..., €, of the n x n identity 
matrix J, as the combination 


Then 


where X;(t) is a particular solution of the simpler equation 
£'(t) = AR(t) + f(DE, f=, €=6). 


An initial trial solution X(t) for ¥’(t) = AX (t)+ f(t)€ can be determined 
from the following initial trial solution rule: 


Assume f(t) is a sum of Euler atoms. Identify independent 
functions whose linear combinations give all derivatives of 
f(t). Let the initial trial solution be a linear combination of 
these functions with undetermined vector coefficients {€ ;}. 


In the well-known scalar case, the trial solution must be modified if its 
terms contain any portion of the general solution to the homogeneous 
equation. In the vector case, if f(t) is a polynomial, then the correc- 
tion rule for the initial trial solution is avoided by assuming the matrix 
A is invertible. This assumption means that r = 0 is not a root of 
det(A — rI) = 0, which prevents the homogenous solution from having 
any polynomial terms. 


The initial vector trial solution is substituted into the differential equa- 


tion to find the undetermined coefficients {¢;}, hence finding a particular 
solution. 


Theorem 29 (Polynomial solutions) 

ket F(t) = Fo Pitt be a polynomial of degree k. Assume A is an n x n 
constant invertible matrix. Then u’ = Ai+/(t)é has a polynomial solution 
uti) = sean Ej hr of degree k with vector coefficients {¢;} given by the 
relations 


OL 


k 
g=— pA Ee, OS5 Sk. 
i=j 
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Theorem 30 (Polynomial x exponential solutions) 

Let g(t) = Die Ps be a polynomial of degree k. Assume A is ann x n 
constant matrix and B = A — al is invertible. Then ti’ = Ati + e%g(t)€ 
has a polynomial-exponential solution u(t) = e@ Lio 55 with vector 
coefficients {¢ ;} given by the relations 


k 
cj=—) piB ec, O<j<k. 
i=j 


Proof of Theorem 29. Substitute u(t) = ae, éj4 into the differential 


equation, then 
2. ee re oe 
j=0 j=0 j=0 
Then terms on the right for 7 = & must add to zero and the others match the 
left side coefficients of t’/j!, giving the relations 


A€y+prE =0, 541 = AG; +p, €. 


Solving these relations recursively gives the formulas 


Ck = —prAeé, 
Ee-1 = — (pr-1A7' + pr A~?) E, 
ee ge ee 


The relations above can be summarized by the formula 
k 
€;=-) pAi* 3g, O<j<k. 
i=j 


The calculation shows that if u(t) = we é;% and €, is given by the last 
formula, then u(t) substituted into the differential equation gives matching 


LHS and RHS. The proof is complete. 


Proof of Theorem 30. Let w(t) = e*’v(t). Then ti’ = Ati +e%g(t)é implies 
v’ = (A-al)¥ + g(t)€. Apply Theorem 29 to v’ = Bv + g(t)€. The proof is 
complete. 
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11.8 Second-order Systems 


A model problem for second order systems is the system of three masses 
coupled by springs studied in section 11.1, equation (6): 


mya (t) = —kyx1(t) + ke[za(t) — 21(8)], 
(1) m2t3(t) = —ke[xa(t) — 21(t)] + kslxa(t) — x(t), 
m3x5(t) = —k3|x3(t) = x2(t)] -_ kaa3(t). 
ky kg keg eg : 
Ms UG Rs UU Figure 22. Three masses 
connected by springs. The masses 
mi, me m3 slide on a frictionless surface. 


In vector-matrix form, this system is a second order system 
ME") = KR) 


where the displacement x, mass matrix MV and stiffness matrix K 
are defined by the formulas 


LY My 0 0 —k, = ko ko 0 
e— x2 , M= 0 m2 0 A ko —ky — kg k3 
x3 0 0 m3 0 kg —kg = ka 


Because M is invertible, the system can always be written as 


x” = AX, A=M'K. 


Converting x” = Ax to u’ = Cu 


Given a second order n x n system X” = AX, define the variable G and 
the 2n x 2n block matrix C’ as follows. 


0 e-(4). e- (448) 


Then each solution X of the second order system X” = AX produces a 
corresponding solution U of the first order system G’ = CU. Similarly, 


each solution i of U’ = CU gives a solution X of ¥” = AX by the 
formula X = diag(J,0)t. 


Euler’s Substitution x = ev 
The fundamental substitution of L. Euler applies to vector-matrix dif- 
ferential systems. In particular, for ¥” = AX, the equation X = ev 


produces the characteristic equation 


det(A — A7I) = 0, 
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and the eigenpair equation 


This eigenpair equation means that (\?,V) is an eigenpair of the matrix 


A. 


Negative eigenvalues of A produce complex conjugate values for 4. 
For instance, \? = —4 implies \ = +27, and then, even though vector V 
has real components, the solution X(t) = eV is a vector with complex 
entries: X(t) = eV = cos(2t)¥ + isin(2t)V. 

Details. Compute xX’ = £ ty = Ney = Ne, Then x = 0° x. I 
x = e¥ is a nonzero solution of ¥” = AX, then \?X = AX holds, 
which is equivalent to A2V = Av. Then (\?,¥) is an eigenpair of A. 
Conversely, if (A?,V) is an eigenpair of A, then the steps reverse to 
obtain \?X = AX, which means that ¥ = e*’¥ is a nonzero solution of 
x= AS, 

By linear algebra, the equation AV = \?¥ has a solution ¥ 4 0 if 
and only if the homogeneous problem (A — A?7J)¥ = 0 has infinitely 
many solutions. Cramer’s Rule implies this event happens exactly when 
det(A — A277) = 0. 


Characteristic Equation for x” = Ax 


The characteristic equation for the n x n second order system ¥” = AX 
will be derived anew from the corresponding 2n x 2n first order system 
u’ = Cu. We will prove the following identity. 


Theorem 31 (Characteristic Equation) 
Let x” = AX be given with n x n constant matrix A. Let i! = Cd be its 
corresponding first order system, where 


fey eSeere 
x! A|0 
Then 
(3) det IAPS) det A= 7h), 


Proof: The method of proof is to verify the product formula 


-M| I EO o| 1 
Al-ar /\OAT[T )) \VA=MF[-ar }° 


Then the determinant product formula applies to give 


(4) det(C — AI) det (art) = act ( poet a 
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Cofactor expansion is applied to give the two identities 
i st ) ~ (—1)" det(A — X72). 


T\0 
aet (Sh r) = det (Spry 


Then (4) implies (3). The proof is complete. 


Solving u’ = Cu and x” = Ax 
Consider the n x n second order system ¥” = AX and its corresponding 


2n x 2n first order system i’ = Cli, where 
: 


0 en($Hh) (3 


Theorem 32 (Eigenanalysis of A and C) 
Let A bea given n xn constant matrix and define the corresponding 2n x 2n 


system by 
a’ = Cu, o= (417), a= ( 


Then 


= Mw, 


Ni Sy 
| 
- 
a 


(6) c-a0 ( )=0 if and only if tee 


Proof: The result is obtained by block multiplication, because 
—AI I 
o-are (244,). 


Theorem 33 (General Solutions of i’ = Cui and xX” = Ax) 
Let A bea given n xn constant matrix and define the corresponding 2n x 2n 


system by 
i’ = Ci, o= (Str). i=(}): 


| ¥2n are independent. 


Assume C has eigenpairs {(Aj,¥;)}521 and yi, 
Let J denote the nxn identity and define w ; = diag(J,0)y¥;, 7 =1,...,2n 
Then u’ = Cti and X” = AX have general solutions 

th = cyer’¥y t+ + cone?"¥ an (2n x 1), 

(t) — ce tw y apres Con O20! W on (n x 1). 


mL Sl 
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Proof: Let ¥;(t) = e*"w;, 7 =1,...,2n. Then X; is a solution of ¥” = AX, 
because X47 (t) = erst (\;)? Wy = AX; (t by Theorem 32. To be verified is the 
independence of the solutions {X pene Let Z; = A;w, and apply Theorem 32 


ee Ww; es a : 
to write y ; = ( Zz. is AW; = NEW Suppose constants aj, ..., Qn are given 
J 


such that bar _, 4X; = 0. Differentiate this relation to give ee ane*Z; =0 


for all t. Set ¢ = 0 in the last summation and combine to obtain See any j = 0. 
Independence of ¥1, ..., Y2n implies that aj = --- = dan = 0. The proof is 
complete. 


Eigenanalysis when A has Negative Eigenvalues. If all cigen- 
values ys of A are negative or zero, then, for some w > 0, eigenvalue pu 
is related to an eigenvalue \ of C by the relation ~ = —w? = \?. Then 

= +wi and w = \/|u]. Consider an eigenpair (—w?, ¥) of the real n x n 
matrix A with w > 0 and let 


ne cy, coswt+cosinwt w>O0, 
~ ) ey +cat w=0. 
Then u(t) = —w?u(t) (both sides are zero for w = 0). It follows 
that X(t) = u(t)V satisfies ¥”(t) = —w?X(t ) and Ax() = u(t)AV = 
—w?X(t). Therefore, ¥(t) = u(t)¥ satisfies ¥”(t) = AX 


Theorem 34 (Eigenanalysis Solution of x” = Ax) 

Let the n x n real matrix A have eigenpairs {(uj,Vj)}7_1. Assume pj = 
—w? with w; > 0, 7 = 1,...,n. Assume that Vj, ..., Vn are linearly 
independent. Then the ener solution of x(t) = Ax(t) is given in terms 


of 2n arbitrary constants aj, ..., Gn, 01, ..., by, by the formula 


i sin wt 
7 t+; J 
(7) a= 3 (% cos w jt + se )s 


in Wt 


This expression uses the limit convention =4. 


Ww 


w=0 


Proof: The text preceding the theorem and superposition establish that X(t) is 
a solution. It only remains to prove that it is the general solution, meaning that 
the arbitrary constants can be assigned to allow any possible initial condition 
X(0) = Xq, ¥’(0) = ¥o. Define the constants uniquely by the relations 


=> n 
xo = dij=1 AgV 5, 
Yo = ye) bj V5, 


which is possible by the assumed independence of the vectors {Vj}_,. Then 
equation (7) implies ¥(0) = )0}_, aj¥j = Xo and X'(0) = Y7_, bj¥ 5 = Yo. 
The proof is complete. 
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11.9 Numerical Methods for Systems 


An initial value problem for a system of two differential equations is 
given by the equations 


'(t 


u(t) 
‘(t) = g(t, x(t), y(¢)), 
(1) < g Y 
y(to) 


A numerical method for (1) is an algorithm that computes an approx- 
imate dot table with first line to, xo, yo. Generally, the dot table has 
equally spaced t-values, two consecutive t-values differing by a constant 
value h # 0, called the step size. To illustrate, if tg = 2, ro = 5, 
yo = 100, then a typical dot table for step size h = 0.1 might look like 


t x y 
2.0 | 5.00 | 100.00 
2.1 | 5.57 | 103.07 
2.2 | 5.62 | 104.10 
2.3 | 5.77 | 102.15 
2.4 | 5.82 | 101.88 
2.5 | 5.96 | 100.55 


Graphics. The dot table represents the data needed to plot a solution 
curve to system (1) in three dimensions (t, x, y) or in two dimensions, 
using a txz-scene or a ty-scene. In all cases, the plot is a simple connect- 
the-dots graphic. 


104 
104 5.8 
100 3D-plot ty-scene ta-scene 
5 aT 58 2 25 


Figure 23. Dot table plots. 
The three dimensional plot is a space curve made directly from the dot ta- 


ble. The ta-scene and the ty-scene are made from the same dot table using 
corresponding data columns. 


Myopic Algorithms. All of the popular algorithms for generating 
a numerical dot table for system (1) are near-sighted, because they 
predict the next line in the dot table from the current dot table line, 
ignoring effects and errors for all other preceding dot table lines. Among 
such algorithms are Euler’s method, Heun’s method and the RK4 
method. 
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Numerical Algorithms: Planar Case 


Stated here without proof are three numerical algorithms for solving two- 
dimensional initial value problems (1). Justification of the formulas is 
obtained from the vector relations in the next subsection. 


Notation. Let to, xo, yo denote the entries of the dot table on a partic- 
ular line. Let h be the increment for the dot table and let top +h, x, y 
stand for the dot table entries on the next line. 
Planar Euler Method. 
= xo + hf (to, £0, Yo), 
y = yo+hg(to, 20, yo): 


Planar Heun Method. 


zy = xo thf (to, 0, yo), 
Y= Yyothg(to, xo, yo), 
zr = «tot+h(f (to, 20, yo) + f(to +h, 21, 491))/2 
y = yot h(g(to, x0, yo) + g(to +h, a1, y1))/2. 


Planar RK4 Method. 


ko o= ae ©0, Yo); 
m, = hg(to, Xo, Yo); 
kg = Ga xo + ki /2, yo + m1/2), 
mz = hg(to +h/2,x0 + k1/2, yo + m1/2), 
kg = hf(to+h/2,x0 + ko/2, yo + m2/2), 
m3 = hg(to + h/2,x9 + k2/2, yo + m2/2), 
kg = hf(to+h,xo + k3, yo + m3), 
m4 = hg(to+h,xo + ks, yo + ms), 
ee no + = (ki + 2k2 + 2k + ka), 

1 
y= yo + & (ma + 2m + 2mz + ma) 


Numerical Algorithms: General Case 


Consider a vector initial value problem 


ii’(t)=F(t,a(t)), (to) = tio. 


Described here are the vector formulas for Euler, Heun and RK4 meth- 
ods. These myopic algorithms predict the next table dot to +h, u from 
the current dot tp, ig. The number of scalar values in a table dot is 
1+, where n is the dimension of the vectors U and F. 
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Vector Euler Method. 
ti = tig + AF (to, tio) 

Vector Heun Method. 
W =o +hF (to, tio), ti = tio + 


(F (to, io) + F (to +h, W)) 


Vector RK4 Method. 


k, = AF (to, tio), 

k, = AF(to+h/2, to +k1/2), 
k, = hE (to +h/2,tio + k2/2), 
k, = hF( lip +k 


a 
a 
j=) 
on 
fon oS 
wl 
bo 
x 
bo 
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